Overlapping Domain Decomposition Methods

81. Overlapping Domain Decomposition Methods#

Let ${Ω_{i}}$ be a non-overlapping decomposition of $Ω$ into open sub-domains, i.e.

\overset{―}{Ω} = ⋃ {\overset{―}{Ω}}_{i}

with

Ω_{i} \cap Ω_{j} = 0 for i \neq j .

Let $H_{i} = diam Ω_{i}$ . To simplify notation we assume that all $H_{i}$ are of the same order of magnitude, and write $H$ .

In addition to $Ω_{i}$ , we define enlarged sub-domains

{\tilde{Ω}}_{i} = U_{c H} (Ω_{i}) = {x \in Ω : dist (x, Ω_{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_{0, D}^{1} (Ω)$ into local spaces supported on the overlapping sub-domains ${\tilde{Ω}}_{i}$ , i.e.

V_{i} = {v \in V_{h} : v = 0 on Ω ∖ {\tilde{Ω}}_{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} with 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^{t h}$ sub-domain ${\tilde{Ω}}_{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} (Ω_{i})$ . The sub-spaces are represented via BitArrays with cleared bits for dofs not belonging to elements outside ${\tilde{Ω}}_{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

σ (C_{D D}^{- 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_{D D} u = inf_{\binom{u = \sum u_{i}}{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 ${ψ_{i}}$ a partition of unity for the domain decomposition ${\tilde{Ω}}_{i}$ iff

support {ψ_{i}} \subset {\tilde{Ω}}_{i}

0 \leq ψ_{i} \leq 1

and

\sum ψ_{i} = 1

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

ψ_{i}^{0} (x) = max {1 - \frac{1}{c H_{i}} dist {x, Ω_{i}}, 0}

These functions are $1$ inside $Ω_{i}$ , $0$ outside ${\tilde{Ω}}_{i}$ , and decay with a maximal derivative

| \nabla ψ_{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:

ψ_{i} (x) := \frac{ψ_{i}^{0} (x)}{\sum ψ_{i}^{0} (x)}

$H^{1}$ -stable quasi-interpolation operators:

We need local interpolation-like operators $ $π_{h} : H^{1} \to 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 π_{h} v ‖_{L_{2} (Ω)} ≺ ‖ \nabla v ‖_{L_{2} (Ω)},

and satisfy the approximation estimate $ $‖ \frac{1}{h} (v - π_{h} v) ‖_{L_{2} (Ω)} ≺ ‖ \nabla v ‖_{L_{2} (Ω)} .$ $

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} = π_{h} (ψ_{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 π_{h} (ψ_{i} u) = π_{h} (\sum ψ_{i} u) = π_{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 π_{h} (ψ_{i} u) ‖_{L_{2}} ≺ ‖ \nabla (ψ_{i} u) ‖_{L_{2}}

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

\begin{array}{rcl} ‖ \nabla u_{i} ‖_{L_{2}} & ≺ & ‖ (\nabla ψ_{i}) u + ψ_{i} \nabla u ‖_{L_{2}} \\ \leq & ‖ \nabla ψ_{i} ‖_{L_{\infty}} ‖ u ‖_{L_{2} ({\tilde{Ω}}_{i})} + ‖ ψ_{i} ‖_{L_{\infty}} ‖ \nabla u ‖_{L_{2} ({\tilde{Ω}}_{i})} \\ ≺ & H^{- 1} ‖ u ‖_{L_{2} ({\tilde{Ω}}_{i})} + ‖ \nabla u ‖_{L_{2} ({\tilde{Ω}}_{i})} \end{array}

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

\begin{array}{rcl} \sum ‖ \nabla u_{i} ‖_{L_{2}}^{2} & ≺ & \sum_{i} {H_{i}^{- 2} ‖ u ‖_{L_{2} ({\tilde{Ω}}_{i})}^{2} + ‖ \nabla u ‖_{L_{2} ({\tilde{Ω}}_{i})}^{2}} \\ ≺ & H^{- 2} ‖ u ‖_{L_{2} (Ω)}^{2} + ‖ \nabla u ‖_{L_{2} (Ω)}^{2} \\ ≺ & H^{- 2} ‖ u ‖_{A}^{2} \end{array}

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_{2 L}$ by the preconditioning action in matrix notation

C_{2 L}^{- 1} = C_{D D}^{- 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 ${\underset{―}{v}}_{H}$ , then the vector $P {\underset{―}{v}}_{H}$ are the coefficients of the same function represented on the fine grid. The matrix $P$ of dimension $\dim V_{h} \times \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.

σ {C_{2 L}^{- 1} A} \subset [c_{1}, c_{2}]

Proof: By the ASM Lemma we have

u^{T} C_{2 L} u = inf_{\binom{u = u_{H} + \sum u_{i}}{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} := π_{H} u

By semi-norm continuity of $π_{H}$ there holds

‖ \nabla u_{H} ‖_{L_{2}} ≺ ‖ \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} ≺ ‖ \nabla u ‖_{L_{2}}^{2}

We proceed as above to decompose $u_{f} = \sum u_{i} = \sum π_{h} (ψ_{i} u_{f})$ , which satisfy the estimates from above

\sum_{i} ‖ \nabla u_{i} ‖^{2} ≺ 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.