Domain Decomposition with minimal overlap

80. Domain Decomposition with minimal overlap#

With domain decomposition (DD) methods one splits a large problem on a big domain $\Omega$ into smaller problems on sub-domains $\Omega_i$. One advantage is that direct factorization on each sub-domain is cheaper than factorization of the global problem. Assume that the global problem is of size $N$, and we use $m$ sub-domains. Assume the (sparse) factorization of a problem of size $n$ has costs $O(n^\alpha)$. Then the advantage is

\[ m \, \big( \tfrac{N}{m} \big)^\alpha = \frac{1}{m^{\alpha-1}} N^\alpha \leq N^\alpha \]

For sparse direct solvers based on nested dissection $\alpha = 1.5$ for 2D problems, and $\alpha = 2$ for 3D.

We sub-divide the unit-square into $m_x \times m_y$ sub-domains. On the lower boundary we assign a boundary condition label:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
mx, my = 3,3
rects = [MoveTo(i/mx,j/mx).Rectangle(1/mx,1/my).Face().bc("mat"+str(i)+str(j)) for i in range(mx) for j in range(my)]
shape = Glue(rects)
shape.edges[Y<0.001].name = "bottom"
Draw (shape);

geo = OCCGeometry(shape,dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
print (mesh.GetMaterials(), mesh.GetBoundaries())
Draw (mesh);

('mat00', 'mat01', 'mat02', 'mat10', 'mat11', 'mat12', 'mat20', 'mat21', 'mat22') ('bottom', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'bottom', 'default', 'default', 'default', 'default', 'default', 'default', 'bottom', 'default', 'default', 'default', 'default', 'default', 'default')

setup a standard problem, Dirichlet b.c. on the bottom boundary:

fes = H1(mesh, order=1, dirichlet='bottom')
print ("ndof =", fes.ndof)
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()

ndof = 488

We build a precontioner as the sum of inverses on sub-domains:

We iterate over the individual sub-domains. The function mesh.Materials returns one volume region, which can be split into a list of regions.
A finite element space can provide a BitArray marking the dofs of a certain region. Then we perform the logical and operation & with the non-Dirichet dofs.
The Inverse(domaindofs) takes the sub-matrix corresponding to the domain dofs, inverts it with a sparse direct solver, and inserts the small inverse into the large by padding with zeros. The costs are $O(N + N_{\text{used}}^\alpha)$, where the constant for the $O(N)$ term is very small.
The final preconditioner is the sum of sub-domain inverses.

pre = None
for dom in mesh.Materials(".*").Split():
    domaindofs = fes.GetDofs(dom) & fes.FreeDofs()
    # print (domaindofs)
    print ("num dofs:", domaindofs.NumSet())
    invi = a.mat.Inverse(domaindofs)
    pre = invi if pre == None else pre + invi
    
from ngsolve.la import EigenValues_Preconditioner
lami = list(EigenValues_Preconditioner(mat=a.mat, pre=pre))
print ("condition number = ", lami[-1]/lami[0])

num dofs: 56
num dofs: 64
num dofs: 64
num dofs: 56
num dofs: 65
num dofs: 65
num dofs: 56
num dofs: 65
num dofs: 65
condition number =  146.02753241632854

To check the solution:

f = LinearForm(x*v*dx).Assemble()
gfu = GridFunction(fes)
from ngsolve.krylovspace import CGSolver
inv = CGSolver(mat=a.mat, pre=pre)
gfu.vec.data = inv * f.vec
print ("Iterations: ", len(inv.errors))
Draw (gfu);

Iterations:  42

Experiment with number of subdomains, and mesh-size. You should observe a condition number $O( h^{-1} H^{-1} )$, where $H$ is the sub-domain size, and $h$ is the element-size.

80.1. Analysis of the method#

We apply the ASM theory, where sub-spaces are

\[ V_i = \operatorname{span} \{ \varphi_i : \text{dof}_i \text{ supported in } \overline{\Omega_i} \} \]

By the ASM Lemma, the preconditioner norm is

\[ \| u \|_{C_{ASM}}^2 = \inf_{u = \sum u_i \atop u_i \in V_i } \sum \| u_i \|_A^2 \]

Since one sub-space $V_i$ has overlap with at most 9 spaces $V_j$, there immediately follows

\[ \| u \|_A^2 \leq 9 \, \| u \|_{C_{ASM}}^2 \]

and

\[ \lambda_{max} (C_{ASM}^{-1} A) \leq 9 \]

For the other bound we have to find a decomposition, and bound the norm:

Lemma Let $u \in H^1(\Omega)$, and $u_i = u$ in $\Omega_i$, and all other dofs set to zero. Then

\[ \| u_i \|_{H_1(\Omega)}^2 \preceq \tfrac{1}{hH} \| u \|_{H^1(\Omega_i)}^2 \]

Proof: Follows from the scaled trace theorem

\[ \| u_i \|_{L_2(\partial \Omega)}^2 \preceq \frac{1}{H} \| u \|_{H^1(\Omega)}^2, \]

and the decay of functions in the first layer of elements outside $\Omega_i$:

\[ \| u_i \|_{H^1(\Omega \setminus \Omega_i)}^2 \preceq \frac{1}{h} \| u_i \|_{L_2(\partial \Omega_i)}^2 \]

$\Box$

Choosing these $u_i$ for every sub-domain does not lead to a decomposition $u = \sum{u_i}$ since dofs at the interface contribute to several sub-space functions. So, we have to divide each coefficient by the number of sub-domains it belongs to, to obtain a decomposition $\tilde u_i$. The estimate is very similar.

gfi = GridFunction(fes)
domi = mesh.Materials(".*").Split()[4]
proj = Projector(fes.GetDofs(domi), True)
gfi.vec.data = proj * gfu.vec
Draw (gfi);

80.2. Adding a coarse grid space#

We build a coarse space of functions being constant inside each sub-domain. Since the functions are continuous, they decay within one layer of elements:

gfcoarse = GridFunction(fes, multidim=len(mesh.Materials(".*").Split()))
mv = gfcoarse.vecs    # a multi-vector
proj = Projector(fes.FreeDofs(),True)
for i,dom in enumerate (mesh.Materials(".*").Split()):
    consti = GridFunction(fes)
    consti.Set(1, definedon=dom)
    proj.Project(consti.vec)
    mv[i].data = consti.vec
Draw (gfcoarse, mesh, animate=True)
print ("Multivector components:",len(mv))

Multivector components: 9

We add the space spanned by these (essentially) piecewise constant functions as additional sub-space, and define a coarse grid correction step:

a0 = InnerProduct(mv, a.mat * mv)
print (a0)
inva0 = BaseMatrix(a0.I)
E0 = BaseMatrix(mv)
coarsegrid = E0 @ inva0 @ E0.T

 23.7946       1       0  2.2975 -0.337032       0       0       0       0
       1 24.9739       1 -0.324899       2 -0.333629       0       0       0
       0       1 16.6868       0 -0.35572       1       0       0       0
  2.2975 -0.324899       0 32.1655       2       0 2.31221 -0.39022       0
 -0.337032       2 -0.35572       2  34.731       2 -0.365027       2 -0.386817
       0 -0.333629       1       0       2 24.7111       0 -0.307914       1
       0       0       0 2.31221 -0.365027       0 23.5794       1       0
       0       0       0 -0.39022       2 -0.307914       1  25.233       1
       0       0       0       0 -0.386817       1       0       1 16.2131

pre2 = pre + coarsegrid
lami = list(EigenValues_Preconditioner(mat=a.mat, pre=pre2))
print ("condition number = ", lami[-1]/lami[0])

condition number =  125.13234921846973

… it does not work very well

We can improve it by scaling the coarse grid functions by the number of domains it belongs to. Thanks to this scaling a globally constant function can be represented within the coarse space.

ones = mv[0].CreateVector()
cnt = mv[0].CreateVector()
cnt[:] = 0
ones[:] = 1
for dom in mesh.Materials(".*").Split():
    cnt += Projector(fes.GetDofs(dom),True)*ones

idiag = DiagonalMatrix(cnt).Inverse(fes.FreeDofs())
for v in mv:
    v.data = (idiag*v).Evaluate()

a0 = InnerProduct(mv, a.mat * mv)
# print (a0)
# print (a0 * Vector( (1,)*len(mv)))
inva0 = BaseMatrix(a0.I)
E0 = BaseMatrix(mv)
coarsegrid = E0 @ inva0 @ E0.T
pre3 = pre + coarsegrid
lami = list(EigenValues_Preconditioner(mat=a.mat, pre=pre3))
print ("condition number = ", lami[-1]/lami[0])

condition number =  22.724988999951798

gfu0 = GridFunction(fes)
gfu0.vec.data = coarsegrid * (a.mat * gfu.vec)
Draw (gfu0);

Exercise: The condition number can be estimated by:

\[ \kappa \{C^{-1} A \} \prec \frac{H}{h} \]

Hint: First, split the global function as $u = u_H + u_f$, where $u_H$ is in domain-wise constant space such that the constant is the mean-value of $u$. Show that

\[ \frac{h}{H} \| \nabla u_f \|_{L_2}^2 + \frac{1}{H^2} \| u_f \|_{L_2}^2 \preceq \| u \|_{H^1}^2 \]

Then, decompose $u_f = \sum u_i$ following the path from above, but using a refined estimates for the trace:

\[ \| u_i \| _{L_2(\partial \Omega)}^2 \preceq \frac{1}{H} \| u_f \|_{L_2(\Omega_i)}^2 + H \| \nabla u_f \|_{L_2(\Omega)}^2 \]

80.3. Graph-based mesh partitioning#

The geometric decomposition of a domain into nice sub-domains as we have done above is usually not feasible. In practice, one often uses graph-based mesh partitioning tools. The mesh connectivity is represented as a graph, where the

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

Draw (mesh);

We build a graph (stored as a list of lists) for the neighbor elements:

nbels = []
for el in mesh.Elements(VOL):
    nbs = []
    for f in el.facets:
        for nb in mesh[f].elements:
            if nb != el:
                nbs.append(nb.nr)
    nbels.append(nbs)
# print (nbels)

Metis is such a graph partitioning tool. Python bindings are provided via the package pymetis. We provide the element neighbor graph as adjacency relation:

import pymetis
ndom = 6
n_cuts, membership = pymetis.part_graph(ndom, adjacency=nbels)

gfdom = GridFunction(L2(mesh,order=0))
gfdom.vec.data = BaseVector(membership)
Draw (gfdom);

A basic method for graph partitioning is spectral bisection: One computes the eigenfunction to the smallest non-zero eigenvalue of the Laplace-operator with Neumann boundary conditions. Two sub-domains are defined for the points where the eigenfunction is positive, or negative. Consider as an example a long rectangle: The eigenfunction is a sine-function in the long direction, and the interface between the sub-domains is as short as possible.

fes = H1(mesh, order=1, dirichlet='bottom')
print ("ndof =", fes.ndof)
print (fes.FreeDofs())
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
gfu = GridFunction(fes)

domdofs = [BitArray(fes.ndof) for i in range(ndom)]
for d in domdofs: d.Clear()
for el in fes.Elements(VOL):
    for d in el.dofs:
        domdofs[membership[el.nr]].Set(d)

ndof = 137
00110000000001111111111111111111111111111111111111
11111111111111111111111111111111111111111111111111
1111111111111111111111111111111111111

pre = None
for dd in domdofs:
    invi = a.mat.Inverse(dd & fes.FreeDofs())
    pre = invi if pre == None else pre + invi
    
from ngsolve.la import EigenValues_Preconditioner
lami = list(EigenValues_Preconditioner(mat=a.mat, pre=pre))
print ("condition number = ", lami[-1]/lami[0])

condition number =  45.10163889987329

80.4. With the coarse space#

gfcoarse = GridFunction(fes, multidim=ndom)
mv = gfcoarse.vecs
proj = Projector(fes.FreeDofs(),True)
for i in range(ndom):
    gfcoarse.Set(IfPos( Norm(gfdom-i)-0.5, 0, 1), mdcomp=i)
    proj.Project(mv[i])

Draw (gfcoarse, animate=True);

a0 = InnerProduct(mv, a.mat * mv)
print (a0)
# print (a0 * Vector( (1,)*len(mv)))
inva0 = BaseMatrix(a0.I)
E0 = BaseMatrix(mv)
coarsegrid = E0 @ inva0 @ E0.T
pre3 = pre + coarsegrid
lami = list(EigenValues_Preconditioner(mat=a.mat, pre=pre3))
print ("condition number =", lami[-1]/lami[0])

 4.46499       0 -2.69243 -1.77256       0       0
       0 11.8079 -2.94996       0 -1.1869       0
 -2.69243 -2.94996 8.13912 -1.73742 -0.759305       0
 -1.77256       0 -1.73742 7.19131 -0.598901 -3.08243
       0 -1.1869 -0.759305 -0.598901  7.8499 -1.7794
       0       0       0 -3.08243 -1.7794 4.86183

condition number = 8.94990000871906