Domain Decomposition with minimal overlap

82. 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:

  1. We iterate over the individual sub-domains. The function mesh.Materials returns one volume region, which can be split into a list of regions.

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

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

  4. 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 =  143.5940776520239

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:  41

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.

82.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);

82.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.4177       1       0 2.30515 -0.343193       0       0       0       0
       1 24.6162       1 -0.32304       2 -0.334919       0       0       0
       0       1 16.3853       0 -0.338679       1       0       0       0
 2.30515 -0.32304       0 31.6869       2       0 2.31213 -0.383486       0
 -0.343193       2 -0.338679       2 34.1546       2 -0.35613       2 -0.375211
       0 -0.334919       1       0       2 24.4771       0 -0.303936       1
       0       0       0 2.31213 -0.35613       0 23.1801       1       0
       0       0       0 -0.383486       2 -0.303936       1 24.6761       1
       0       0       0       0 -0.375211       1       0       1 15.9986

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

condition number =  123.20090512279653

… 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.309944887784315

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 \]

82.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)

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[16], line 1
----> 1 import pymetis
      2 ndom = 6
      3 n_cuts, membership = pymetis.part_graph(ndom, adjacency=nbels)

ModuleNotFoundError: No module named 'pymetis'

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

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[17], line 2
      1 gfdom = GridFunction(L2(mesh,order=0))
----> 2 gfdom.vec.data = BaseVector(membership)
      3 Draw (gfdom);

NameError: name 'membership' is not defined

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 = 136
0: 00110000000001111111111111111111111111111111111111
50: 11111111111111111111111111111111111111111111111111
100: 111111111111111111111111111111111111

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[18], line 8
      4 u,v = fes.TnT()
      5 a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
      6 gfu = GridFunction(fes)
      7 
----> 8 domdofs = [BitArray(fes.ndof) for i in range(ndom)]
      9 for d in domdofs: d.Clear()
     10 for el in fes.Elements(VOL):
     11     for d in el.dofs:

NameError: name 'ndom' is not defined

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])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[19], line 2
      1 pre = None
----> 2 for dd in domdofs:
      3     invi = a.mat.Inverse(dd & fes.FreeDofs())
      4     pre = invi if pre == None else pre + invi
      5 

NameError: name 'domdofs' is not defined

82.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);

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[20], line 1
----> 1 gfcoarse = GridFunction(fes, multidim=ndom)
      2 mv = gfcoarse.vecs
      3 proj = Projector(fes.FreeDofs(),True)
      4 for i in range(ndom):

NameError: name 'ndom' is not defined

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])

---------------------------------------------------------------------------
NgException                               Traceback (most recent call last)
Cell In[21], line 1
----> 1 a0 = InnerProduct(mv, a.mat * mv)
      2 print (a0)
      3 # print (a0 * Vector( (1,)*len(mv)))
      4 inva0 = BaseMatrix(a0.I)

NgException: BaseVector::InnerProduct: size of me = 488 != size of other = 136
Contents