BDDC - Preconditioner

99. BDDC - Preconditioner#

The idea of the BDDC preconditioner grew out of FETI-DP methods. It was invented in Clark Dormann: A preconditioner for substructuring based on constrained energy minimization, 2003. BDDC stands for Balancing Domain Decomposition with Constraints.

It is a preconditioner for non-overlapping domain decomposition. Unlike the FETI method, it can be used to solve the system for the primal variable only, and allows to use conjugate gradient iteration for the positive definite system.

The idea is to define a problem on a larger space, where the matrix is cheaper to invert. Solve the problem on the larger space, and project down the (artificial) solution back to the original space. Before, the original right hand side was lifted into the larger space by the adjoint of the projection.

Theorem: Let

$A(.,.) : V \times V \rightarrow {\mathbb R}$ be an symmetric and elliptic bilinear-form
$\widetilde V \supset V$ a larger space
$\widetilde A(.,.) : \widetilde V \times \widetilde V \rightarrow {\mathbb R}$ be an extension of $A(.,.)$
Let $R : \widetilde V \rightarrow V$ be a projection

The fictitious space preconditioner is $$ C^{-1} := R \widetilde A^{-1} R^T $$

Its spectrum can can be estimated as

\[ \sigma (C^{-1} A) \subset [1, \| R \|_{(\widetilde V, \| \cdot \|_\widetilde A) \rightarrow (V, \| \cdot \|_A)}^2] \]

Proof: By the additive Schwarz lemma the preconditioner $C$ can be represented as:

\[ \| u \|_C^2 = \inf_{\tilde u \in \widetilde V \atop R \tilde u = u} \; \| \tilde u \|_{\widetilde A}^2 \]

Since $V \subset \widetilde V$, and $\widetilde A(.,.)$ coincides with $A(.,.)$ on $V$, we immediately have:

\[ \| u \|_C^2 \leq \| u \|_A^2 \]

The other side follows from

\[\begin{eqnarray*} \sup_u \frac{\| u \|_A^2 }{ \| u \|_C^2 } & = & \sup_u \frac{\| u \|_A^2 }{ \inf_{\tilde u : R \tilde u = u} \| \tilde u \|_{\tilde A^2 }} \\ & = & \sup_u \sup_{\tilde u : R \tilde u = u} \frac{ \| u \|_A^2 }{ \| \tilde u \|_{\tilde A}^2} \\ & = & \sup_u \sup_{\tilde u : R \tilde u = u} \frac{ \| R \tilde u \|_A^2 }{ \| \tilde u \|_{\tilde A}^2} \\ & = & \sup_{\tilde u \in widetilde V} \frac{ \| R \tilde u \|_A^2 }{ \| \tilde u \|_{\tilde A}^2} \\ & = & \| R \|_{(\widetilde V, \| \cdot \|_\widetilde A) \rightarrow (V, \| \cdot \|_A)}^2 \end{eqnarray*}\]

We apply the theorem as follows:

the original space is the globally continuous finite element space
the larger space consists of sub-domain wise continuous finite element spaces, which are globally continuous only in the sub-domain vertices (scroll down to see such functions)
The global system is cheaper to invert. One can eliminate all variables except at the vertices locally. The resulting Schur complement (for the vertices) form a global system witch is much smaller than the original system. This global vertex-system brings the necessary global interaction to the method

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *
shape = None
mx, my = 3,3
for i in range(mx): 
    for j in range(my):
        rect = MoveTo(i/mx, j/mx).Rectangle(1/mx, 1/my).Face()
        rect.faces.name='mat'+str(i)+str(j)
        shape = Glue([shape,rect]) if shape else rect

shape.edges[Y<1e-4].name = "bot"
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.04))
Draw (mesh);

fes = None
for dom in mesh.Materials('.*').Split():
    fesi = Compress(H1(mesh, definedon=dom, dirichlet="bot"))
    fes = fes * fesi if fes else fesi

fesVertex = H1(mesh, definedon=mesh.BBoundaries('.*'))
fes = fes * fesVertex

freedofs = fes.FreeDofs()
for el in fes.Elements(BBND):
    freedofs.Set(el.dofs[-1])

u, v = fes.TnT()
domtrial = {} 
domtest = {}
for nr,dom in enumerate (mesh.Materials('.*').Split()):
    domtrial[dom] = u[nr]
    domtest[dom] = v[nr]

a = BilinearForm(fes)
f = LinearForm(fes)

uvert,vvert = u[-1], v[-1]
import ngsolve.comp
dVert = ngsolve.comp.DifferentialSymbol(BBND)

for dom in mesh.Materials('.*').Split():
    ui, vi = domtrial[dom], domtest[dom]
    a += grad(ui)*grad(vi)*dx 
    a += 1e6*(ui-uvert)*(vi-vvert)*dVert(dom.Boundaries().Boundaries())
    f += y*x*vi*dx
            
a.Assemble()
f.Assemble()

gfu = GridFunction(fes)

gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec

gftot = CF ( list(gfu.components) )
Draw(gftot, mesh);