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:

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

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

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)

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

 1.00229
 1.04141
 1.13674
 ...
  7.35531
 7.49161
 7.66495

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

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

 1.00137
 1.03997
 1.11641
 ...
  6.92606
 7.01828
 7.29664

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.

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

 1.00412
 1.09465
   1.334
 ...
  31.0839
 33.0367
 33.6481

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:

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)

18.170358091446623

This is a slight improvement from the default.

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

27.352602492150023

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.

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 =  1154741

WARNING: kwarg 'coarsetype' is an undocumented flags option for class <class 'ngsolve.comp.Preconditioner'>, maybe there is a typo?

H1AMG: level = 0, num_edges = 57801, nv = 1154741
H1AMG: level = 1, num_edges = 30592, nv = 4795
H1AMG: level = 2, num_edges = 16407, nv = 2509
H1AMG: level = 3, num_edges = 8552, nv = 1315
H1AMG: level = 4, num_edges = 4452, nv = 692
H1AMG: level = 5, num_edges = 2265, nv = 361
H1AMG: level = 6, num_edges = 1140, nv = 192
H1AMG: level = 7, num_edges = 574, nv = 102
H1AMG: level = 8, num_edges = 270, nv = 54
H1AMG: level = 9, num_edges = 121, nv = 28
H1AMG: level = 10, num_edges = 57, nv = 15
CG iteration 1, residual = 0.7043503007220652
CG iteration 2, residual = 0.2806108630724647
CG iteration 3, residual = 0.2571230982155469
CG iteration 4, residual = 0.2978589697772882
CG iteration 5, residual = 0.26029187831704664
CG iteration 6, residual = 0.20505150436287917
CG iteration 7, residual = 0.16198117954140395
CG iteration 8, residual = 0.1361015811683436
CG iteration 9, residual = 0.10123085436239918
CG iteration 10, residual = 0.068624149488305
CG iteration 11, residual = 0.05096951764653183
CG iteration 12, residual = 0.041118407847676196
CG iteration 13, residual = 0.03398634520396559
CG iteration 14, residual = 0.026760005142223934
CG iteration 15, residual = 0.018926794946863425
CG iteration 16, residual = 0.014671034379641075
CG iteration 17, residual = 0.01209257890350033
CG iteration 18, residual = 0.009100371585157615
CG iteration 19, residual = 0.006812240674031471
CG iteration 20, residual = 0.004913766208512962
CG iteration 21, residual = 0.003919240361008209
CG iteration 22, residual = 0.0030023276232985837
CG iteration 23, residual = 0.0022813322987511697
CG iteration 24, residual = 0.0017029999959415818
CG iteration 25, residual = 0.0013471279759169437
CG iteration 26, residual = 0.0010667207279254156
CG iteration 27, residual = 0.0008395212065148367
CG iteration 28, residual = 0.0006316755388500332
CG iteration 29, residual = 0.00048795188939796297
CG iteration 30, residual = 0.00037782005830912295
CG iteration 31, residual = 0.00029054307494693515
CG iteration 32, residual = 0.00023026314035244985
CG iteration 33, residual = 0.00017807239766605973
CG iteration 34, residual = 0.00013613775587554015
CG iteration 35, residual = 0.00010788658795989151
CG iteration 36, residual = 8.441304834933107e-05
CG iteration 37, residual = 6.529501167682916e-05
CG iteration 38, residual = 4.994735406139828e-05
CG iteration 39, residual = 3.8493987960052186e-05
CG iteration 40, residual = 3.062063636692641e-05
CG iteration 41, residual = 2.367809919415635e-05
CG iteration 42, residual = 1.8299137201962535e-05
CG iteration 43, residual = 1.4349281485019734e-05
CG iteration 44, residual = 1.0973177553334421e-05
CG iteration 45, residual = 8.931950138092949e-06
CG iteration 46, residual = 7.88208233378162e-06
CG iteration 47, residual = 5.818669864540491e-06
CG iteration 48, residual = 4.4206613691773024e-06
CG iteration 49, residual = 3.452922364422291e-06
CG iteration 50, residual = 2.7440803357712447e-06
CG iteration 51, residual = 2.072639611652896e-06
CG iteration 52, residual = 1.5739217780609395e-06
CG iteration 53, residual = 1.1958288352524757e-06
CG iteration 54, residual = 9.250110962831961e-07
CG iteration 55, residual = 7.309849270710021e-07
CG iteration 56, residual = 5.812597025686664e-07
CG iteration 57, residual = 4.378365219084341e-07
CG iteration 58, residual = 3.3964175579654306e-07
CG iteration 59, residual = 2.5894784528779616e-07
CG iteration 60, residual = 2.01668968614017e-07
CG iteration 61, residual = 1.6062498075548296e-07
CG iteration 62, residual = 1.2615730769350752e-07
CG iteration 63, residual = 9.644037778046905e-08
CG iteration 64, residual = 7.386281083587735e-08
CG iteration 65, residual = 5.67968105987631e-08
CG iteration 66, residual = 4.3560976622339744e-08
CG iteration 67, residual = 3.427920837400871e-08
CG iteration 68, residual = 2.632550278258314e-08
CG iteration 69, residual = 2.074146447389145e-08
CG iteration 70, residual = 1.6311204230090576e-08
CG iteration 71, residual = 1.2729626873837223e-08
CG iteration 72, residual = 9.671932986101752e-09
CG iteration 73, residual = 7.584899636291768e-09
CG iteration 74, residual = 5.88734386185073e-09
CG iteration 75, residual = 4.5098375419924126e-09
CG iteration 76, residual = 3.4675859002578044e-09
CG iteration 77, residual = 2.6711371724206638e-09
CG iteration 78, residual = 2.0776233426675343e-09
CG iteration 79, residual = 1.634584027013332e-09
CG iteration 80, residual = 1.2864300336585845e-09
CG iteration 81, residual = 9.713822015985417e-10
CG iteration 82, residual = 7.682108493625249e-10
CG iteration 83, residual = 6.418212198694404e-10
CG iteration 84, residual = 5.509418093417826e-10
CG iteration 85, residual = 4.0073391807306437e-10
CG iteration 86, residual = 2.9800211301491927e-10
CG iteration 87, residual = 2.261095079735559e-10
CG iteration 88, residual = 1.742677600595856e-10
CG iteration 89, residual = 1.3526528574880763e-10
CG iteration 90, residual = 1.0715192061577473e-10
CG iteration 91, residual = 8.502477796398921e-11
CG iteration 92, residual = 6.339968080693091e-11
CG iteration 93, residual = 4.836023265334509e-11
CG iteration 94, residual = 3.6868423594093425e-11
CG iteration 95, residual = 2.789384134808953e-11
CG iteration 96, residual = 2.1565843842843582e-11
CG iteration 97, residual = 1.642808903219412e-11
CG iteration 98, residual = 1.2498302127522263e-11
CG iteration 99, residual = 9.742134105010631e-12
CG iteration 100, residual = 7.862021306246275e-12
CG iteration 101, residual = 6.056733810186065e-12
CG iteration 102, residual = 4.65990318386506e-12
CG iteration 103, residual = 3.5072706075568113e-12
CG iteration 104, residual = 2.6759263128335983e-12
CG iteration 105, residual = 2.040890461630778e-12
CG iteration 106, residual = 1.5638846641457056e-12
CG iteration 107, residual = 1.1985360285873696e-12
CG iteration 108, residual = 9.275570368120818e-13
CG iteration 109, residual = 7.37642471665576e-13
CG iteration 110, residual = 5.790744649109994e-13

BaseWebGuiScene

Postscript

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

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 = 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],deformation=True)

Newton iteration  0
err =  0.39346141669994983
Newton iteration  1
err =  1.8488788182878563e-15

BaseWebGuiScene