# 2.1.3 Element-wise BDDC Preconditioner¶

The element-wise BDDC 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.

A simple description of the element-wise BDDC (Balancing Domain Decomposition preconditioner with Constraints) can be given in the context of Lagrange finite element space $$V$$: 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. NGSolve classifies dofs into (condensable) LOCAL_DOFs and a remainder that consists of these types:

WIREBASKET_DOF
INTERFACE_DOF

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

In [1]:

import netgen.gui
# %gui tk
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)

In [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.

In [3]:

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

In [4]:

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

 1.00145
1.05935
1.23239
...
21.2736
22.7821
22.9382



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

In [5]:

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

 1.00026
1.03169
1.09454
...
4.10312
4.13286
4.24601



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). This is the default in 3D.

The complete situation is a bit more complex due to the fact these options can take three values (True, False, Undefined), can be combined, and space dimension can be 2 or 3. So 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.

In [6]:

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

 1.00121
1.0817
1.23231
...
25.7229
25.723
27.2163



Clearly, the conditioning became less favorable compared to the default.

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

You can view the types of dofs in a finite element space as follows.

In [7]:

fes = H1(mesh, order=3)
print ("number vertices =", mesh.nv)
print ("dofs of first edge: ", fes.GetDofNrs(NodeId(EDGE,0)))
for i in range(fes.ndof):
print ("ct[",i,"] = ", fes.CouplingType(i) )

number vertices = 21
dofs of first edge:  (21, 22)
ct[ 22 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 24 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 26 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 28 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 30 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 32 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 34 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 36 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 38 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 40 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 42 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 44 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 46 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 48 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 50 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 52 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 54 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 56 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 58 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 60 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 62 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 64 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 66 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 68 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 70 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 72 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 74 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 76 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 78 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 80 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 82 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 84 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 86 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 88 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 90 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 92 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 94 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 96 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 98 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 100 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 102 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 104 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 106 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 108 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 110 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 112 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 114 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 116 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 118 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 120 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 122 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 124 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 126 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 128 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 130 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 132 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 134 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 136 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 138 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 140 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 142 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 144 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 146 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 148 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 150 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 152 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 153 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 154 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 155 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 156 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 157 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 158 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 159 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 160 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 161 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 162 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 163 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 164 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 165 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 166 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 167 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 168 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 169 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 170 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 171 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 172 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 173 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 174 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 175 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 176 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 177 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 178 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 179 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 180 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 181 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 182 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 183 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 184 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 185 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 186 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 187 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 188 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 189 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 190 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 191 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 192 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 193 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 194 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 195 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 196 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 197 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 198 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 199 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 200 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 201 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 202 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 203 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 204 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 205 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 206 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 207 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 208 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 209 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 210 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 211 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 212 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 213 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 214 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 215 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 216 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 217 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 218 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 219 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 220 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 221 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 222 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 223 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 224 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 225 ] =  COUPLING_TYPE.INTERFACE_DOF
ct[ 226 ] =  COUPLING_TYPE.INTERFACE_DOF

In [8]:

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

Out[8]:

9.876864325521144


We may modify the dof types as follows.

In [9]:

fes = H1(mesh, order=8)
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

a = BilinearForm(fes, eliminate_internal=True)
c = Preconditioner(a, "bddc")
a.Assemble()

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

Out[9]:

6.717426419118481

In [ ]: