Overlapping Domain Decomposition Methods

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

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

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.06687894741407885 2.9999999999993516
kappa= 44.85716531130414

We observe the following behaviour:

the condition number is independent of the refinement level
the condition number grows with the number of sub-domains

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

67.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()

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.0427258429588364 3.692782966558739
kappa= 3.5414706478167903

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

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