81. Overlapping Domain Decomposition Methods

Let \(\{ \Omega_i \}\) be a non-overlapping decomposition of \(\Omega\) into open sub-domains, i.e.

\[ \overline \Omega = \bigcup \ \overline \Omega_i \]

with

\[ \Omega_i \cap \Omega_j = 0 \quad \text{ for } i \neq j. \]

Let \(H_i = \operatorname{diam} \Omega_i\). To simplify notation we assume that all \(H_i\) are of the same order of magnitude, and write \(H\).

In addition to \(\Omega_i\), we define enlarged sub-domains

\[ \widetilde \Omega_i = {U}_{cH} (\Omega_i) = \{ x \in \Omega : \operatorname{dist} (x,\Omega_i) < c H \}, \]

where the constant \(c\) is of order one.

The overlapping domain decomposition method is a sub-space decomposition of the global finite element space \(V_h \subset H^1_{0,D}(\Omega)\) into local spaces supported on the overlapping sub-domains \(\widetilde \Omega_i\), i.e.

\[ V_i = \{ v \in V_h : v = 0 \text{ on } \Omega \setminus \widetilde \Omega_i \}. \]

The overlapping domain decomposition method is an additive Schwarz method with the space splitting

\[ V_h = \sum V_i \]

We are solving Dirichlet problems on the overlapping sub-domains, and adding up the solutions: The preconditioning action

\[ C^{-1} : r \mapsto w \]

is defined as

\[ w = \sum w_i \quad \text{with} \quad w_i \in V_i : A(w_i, v_i) = r(v_i) \; \forall \, v_i \in V_i \]

81.1. Experiments with overlapping DD

from ngsolve import *
from ngsolve.webgui import Draw

We start from a coarse grid, where we define the overlapping sub-domains as vertex patches: All elements connected with vertex i form the \(i^{th}\) sub-domain \(\widetilde \Omega_i\). To visualize it, we define an \(L_2\) finite element function with order=0.

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fesdom = L2(mesh, order=0, autoupdate=True)
fes = H1(mesh, order=1, dirichlet="bottom", autoupdate=True)
gfdom = GridFunction(fesdom, multidim=mesh.nv, nested=True, autoupdate=True)
for el in mesh.Elements(VOL):
    for v in el.vertices:
        gfdom.vecs[v.nr][el.nr] = 1
        
for l in range(2):
    mesh.Refine()

Draw (gfdom, animate=True);

copy mesh complete

Setup the problem. The sub-spaces are now finite element functions which are in \(H_0^1(\Omega_i)\). The sub-spaces are represented via BitArrays with cleared bits for dofs not belonging to elements outside \(\widetilde \Omega_i\):

u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()

pre = None
for v in gfdom.vecs:
    mask = BitArray(fes.FreeDofs())
    for el in fes.Elements(VOL):
        if v[el.nr] == 0:
            for d in el.dofs:
                mask[d] = False
    # print (mask)
    invi = a.mat.Inverse(freedofs=mask, inverse="sparsecholesky")
    pre = pre+invi if pre else invi     

from ngsolve.la import EigenValues_Preconditioner
lam = list(EigenValues_Preconditioner(a.mat, pre))
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 0.06687894741407863 2.99999999999935
kappa= 44.857165311304264

We observe the following behaviour:

• the condition number is independent of the refinement level

• the condition number grows with the number of sub-domains

81.2. Analysis of the DD preconditioner

Theorem: The eigenvalues of the domain decomposition preconditioner are bounded as

\[ \sigma(C_{DD}^{-1} A) \subset [c_1 H^2, c_2 ] \]

Outline of the proof: We apply the Additive Schwarz Lemma. The ASM lemma provides the representation

\[ u^T C_{DD} u = \inf_{u = \sum u_i \atop u_i \in V_i} \sum_i \| u_i \|_A^2 \]

By assuming shape regularity of the domains, only a small, fixed number of sub-domains overlap. From that we get the constant upper bound \(c_2\) for the spectrum. For the lower bound \(c_1 H^2\) we have to define explicitly a decomposition of a finite element function \(u\) into sub-space finite element functions. This requires some technical tools.

Partition on unity:

We call the set of functions \(\{ \psi_i \}\) a partition of unity for the domain decomposition \(\widetilde \Omega_i\) iff

\[ \text{support} \{ \psi_i \} \subset \widetilde \Omega_i \]

\[ 0 \leq \psi_i \leq 1 \]

and

\[ \sum \psi_i = 1 \]

For the definition of the overlapping domain decomposition we can define a partition of unity as follows:

\[ \psi_i^0(x) = \max \{1 - \frac{1}{cH_i} \operatorname{dist} \{ x, \Omega_i \}, 0 \} \]

These functions are \(1\) inside \(\Omega_i\), \(0\) outside \(\widetilde \Omega_i\), and decay with a maximal derivative

\[ | \nabla \psi_i^0 |_{\infty} \leq \frac{1}{c H_i}. \]

However, these function do not sum up to 1, so we rescale by their sum to obtain a partition of unity:

\[ \psi_i(x) := \frac{\psi_i^0(x)} {\sum \psi_i^0(x) } \]

\(H^1\)-stable quasi-interpolation operators:

We need local interpolation-like operators $\( \pi_h : H^1 \rightarrow V_h \)$

which are projectors onto \(V_h\), i.e. \(I_h v_h = v_h\) for \(v_h \in V_h\), are continuous in the \(H^1\) semi-norm

\[ \| \nabla \pi_h v \|_{L_2(\Omega)} \prec \| \nabla v \|_{L_2(\Omega)}, \]

and satisfy the approximation estimate $\( \| \tfrac{1}{h} (v - \pi_h v) \|_{L_2(\Omega)} \prec \| \nabla v \|_{L_2(\Omega)}. \)$

We cannot use simple nodal interpolation since this is not well-defined for \(H^1\) functions (for two or more space dimensions), but we can use a Clément-type quasi-interpolation operator: Such an operator defines nodal values by some local averaging.

We have now the tools to conduct the proof for the sub-space decomposition:

For a given \(u \in V_h\) we define

\[ u_i = \pi_h (\psi_i u) \]

These \(u_i\) are in \(V_i\), i.e. are localized finite element functions.

By linearity of the quasi-interpolation, the partition of unity, and the projection property they are a decomposition of the given function \(u\):

\[ \sum u_i = \sum \pi_h (\psi_i u) = \pi_h \big( \sum \psi_i u \big) = \pi_h u = u \]

Finally, we have to bound the norm of the decomposition. First, we use the boundedness of the quasi-interpolation in the \(H^1\) semi-norm:

\[ \| \nabla u_i \|_{L_2} = \| \nabla \pi_h (\psi_i u) \|_{L_2} \prec \| \nabla (\psi_i u) \|_{L_2} \]

We continue with the product rule, and bounds for the pu-functions and their derivatives:

\[\begin{eqnarray*} \| \nabla u_i \|_{L_2} & \prec & \| (\nabla \psi_i) u + \psi_i \nabla u \|_{L_2} \\ & \leq & \| \nabla \psi_i \|_{L_\infty} \| u \|_{L_2(\widetilde \Omega_i)} + \| \psi_i \|_{L_\infty} \| \nabla u \|_{L_2(\widetilde \Omega_i)} \\ & \prec & H^{-1} \, \| u \|_{L_2(\widetilde \Omega_i)} + \| \nabla u \|_{L_2(\widetilde \Omega_i)} \end{eqnarray*}\]

Finally, we sum over all sub-domains. We use that only a small, fixed number of sub-domains overlap:

\[\begin{eqnarray*} \sum \| \nabla u_i \|_{L_2}^2 & \prec & \sum_i \Big\{ H_i^{-2} \| u \|_{L_2(\widetilde \Omega_i)}^2 + \| \nabla u \|_{L_2(\widetilde \Omega_i)}^2 \Big\} \\ & \prec & H^{-2} \| u \|_{L_2(\Omega)}^2 + \| \nabla u \|_{L_2(\Omega)}^2 \\ & \prec & H^{-2} \| u \|_A^2 \end{eqnarray*}\]

In the last step we used Friedrichs’ inequality to bound the \(L_2\)-norm by the \(H^1\) semi-norm. The factor \(H^{-2}\) we paid for the derivative of the pu-functions shows up in the final result.

81.3. Overlapping DD Methods with coarse grid

Now we add an additional coarse grid correction step. We assume that we have an additional finite element space \(V_H\), where the element-size is in the order of the sub-domain size \(H\). We assume that \(V_H \subset V_h\).

We define this two-level preconditioner \(C_{2L}\) by the preconditioning action in matrix notation

\[ C_{2L}^{-1} = C_{DD}^{-1} + P A_H^{-1} P^T, \]

where \(A_H\) is the discretization matrix on the coarse grid, and \(P\) is the so called prolongation matrix. If \(v_H\) is some finite element function in the coarse space \(V_H\) with coefficient vector \(\underline v_H\), then the vector \(P \underline v_H\) are the coefficients of the same function represented on the fine grid. The matrix \(P\) of dimension \(\operatorname{dim} V_h \times \operatorname{dim} V_H\) will not be stored as a dense matrix, but the matrix vector product and its transpose are implemented as linear operators. They are called prolongation and restriction operators.

from ngsolve import *
from ngsolve.webgui import Draw

Same as before, but now we assemble the bilinear-form already before mesh refinement, and invert the (small) coarse grid matrix by a direct solver:

mesh = Mesh(unit_square.GenerateMesh(maxh=0.3))
fesdom = L2(mesh, order=0, autoupdate=True)
fes = H1(mesh, order=1, dirichlet="bottom", autoupdate=True)
gfdom = GridFunction(fesdom, multidim=mesh.nv, nested=True, autoupdate=True)
for el in mesh.Elements(VOL):
    for v in el.vertices:
        gfdom.vecs[v.nr][el.nr] = 1
 
u,v = fes.TnT()
a = BilinearForm(grad(u)*grad(v)*dx).Assemble()
a0inv = a.mat.Inverse(fes.FreeDofs())

levels = 2
for l in range(levels):
    mesh.Refine()

copy mesh complete

the local inverses:

a.Assemble()

pre = None
for v in gfdom.vecs:
    mask = BitArray(fes.FreeDofs())
    for el in fes.Elements(VOL):
        if v[el.nr] == 0:
            for d in el.dofs:
                mask[d] = False

    invi = a.mat.Inverse(freedofs=mask, inverse="sparsecholesky")
    pre = pre+invi if pre else invi     

We have the prolongation as a function call, prolongating a vector in-place. Prolongate(l) prolongates a vector from level l-1 to level l. The transpose of the prolongation matrix is the restriction matrix, again implemented as an in-place operation:

class ProlongationOp(BaseMatrix):
    def __init__ (self, fes, levels, nc):
        super(ProlongationOp, self).__init__()
        self.fes = fes
        self.levels = levels
        self.nc = nc
        
    def Mult (self, x, y):
        y[IntRange(0, self.nc)] = x
        for l in range(self.levels):
            fes.Prolongation().Prolongate(l+1, y)

    def MultTrans (self, y, x):
        hy = y.CreateVector()
        hy.data = y
        for l in range(self.levels, 0, -1):
            fes.Prolongation().Restrict(l, hy)
        x.data = hy[IntRange(0, self.nc)]

    def Height (self):
        return self.fes.ndof
    def Width (self):
        return self.fes.ndof

P = ProlongationOp (fes, levels, a0inv.height)

C2l = P @ a0inv @ P.T + pre

lam = list(EigenValues_Preconditioner(a.mat, C2l))
print ("lammin, lammax=", lam[0], lam[-1])
print ("kappa=", lam[-1]/lam[0])

lammin, lammax= 1.0427258429588369 3.692782966558738
kappa= 3.5414706478167877

81.4. Analysis of the 2-level method:

We apply now the ASM - theory to the sub-space splitting

\[ V_h = V_H + \sum V_i \]

Theorem: The condition number of the 2-level method is bounded, i.e.

\[ \sigma \{ C_{2L}^{-1} A \} \subset [c_1, c_2] \]

Proof: By the ASM Lemma we have

\[ u^T C_{2L} u = \inf_{u = u_H + \sum u_i \atop u_H \in V_H, u_i \in V_i} \, \| u_H \|_A^2 + \sum_i \| u_i \|_A^2. \]

Again, only a fixed number of sub-spaces overlap, which immediately proves the constant \(c_2\).

We have to come up with an explicit decomposition \(u = u_H + \sum u_i\). First we fix the coarse grid function \(u_H\) by Clément interpolation onto the coarse-grid:

\[ u_H := \pi_H u \]

By semi-norm continuity of \(\pi_H\) there holds

\[ \| \nabla u_H \|_{L_2} \prec \| \nabla u \|_{L_2}, \]

and for the rest \(u_f = u - u_H\) there holds \(H^1\)-stability as well as good approximation in \(L_2\)-norm:

\[ \| \nabla u_f \|_{L_2}^2 + \| H^{-1} u_f \|_{L_2}^2 \prec \| \nabla u \|_{L_2}^2 \]

We proceed as above to decompose \(u_f = \sum u_i = \sum \pi_h (\psi_i u_f)\), which satisfy the estimates from above

\[ \sum_i \| \nabla u_i \|^2 \prec H^{-2} \| u_f \|_{L_2}^2 + \| \nabla u_f \|_{L_2}^2 \]

The inverse factor \(H^{-2}\) we pay due to the partition of unity, we get back from the approximation error estimate of the coarse grid interpolant.

81.5. Comparison to DD with minimal overlap

In comparison to the small overlap in the previous section, we have now an overlap of a fixed width, independent of the mesh-size. It is more expensive, but the condition number does not deteriorate for smaller mesh-size. The second difference is that we have now a coarse space containing the continuous, piecewise linear finite element functions. Thanks to this coarse space, we have a good approximation property of coarse grid functions.