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.