On the Scalability of the Schwarz Method

TL;DR

This paper investigates the convergence and scalability of the one-level Parallel Schwarz method for complex domain decomposition problems, especially those with interior subdomains and multiple overlaps, extending existing theoretical results.

Contribution

It develops a systematic framework to analyze the Schwarz method's convergence for interior subdomains with arbitrary overlaps, broadening previous limited scenarios.

Findings

Bounded the norm of the Schwarz iteration operator for interior subdomains.

Extended scalability analysis to multiple overlapping subdomains.

Provided theoretical insights applicable to solvation models in chemistry.

Abstract

In this article, we analyse the convergence behaviour and scalability properties of the one-level Parallel Schwarz method (PSM) for domain decomposition problems in which the boundaries of many subdomains lie in the interior of the global domain. Such problems arise, for instance, in solvation models in computational chemistry. Existing results on the scalability of the one-level PSM are limited to situations where each subdomain has access to the external boundary, and at most only two subdomains have a common overlap. We develop a systematic framework that allows us to bound the norm of the Schwarz iteration operator for domain decomposition problems in which subdomains may be completely embedded in the interior of the global domain and an arbitrary number of subdomains may have a common overlap.

Figures

Equations

$$\begin{array}{rl}\Delta u& =0\text{in }\Omega ,\\ u& =g\text{on }\partial \Omega .\end{array}$$

$${\Gamma }_i\cap {\Omega }_k$$

$${\Gamma }_j\cap {\Omega }_k$$

$${\Gamma }_j^{\text{ext}}:=\partial {\Omega }_j\cap \partial \Omega ,{\Gamma }_j^{\mathrm{int}}:=\overline{\partial {\Omega }_j\setminus {\Gamma }_j^{\text{ext}}}.$$

$$N_j:=\{i\in \mathbb{N},i\ne j:{\Omega }_i\cap {\Omega }_j\ne \varnothing \}.$$

$$N_{jk}:=(N_j\cap N_k)\cup \{0\}.$$

$$\Gamma_{j}^{k,i}:=\begin{cases}\text{int}\left\{x\in\partial\Omega_{j}\colon x\in\Omega_{k}\cap\Omega_{i}\right\}\quad&\text{ if }i\neq 0,\\ \text{int}\left\{x\in\partial\Omega_{j}\colon(x\in\Omega_{k})\land\big{(}x\notin\Omega_{\ell}\leavevmode\nobreak\ \forall\ell\in N_{j}\text{ such that }\ell\neq k\big{)}\right\}\quad&\text{ if }i=0.\end{cases}$$

$${\Gamma }_j^{\mathrm{int}}=\overline{{\cup }_{k\in N_j}{\cup }_{i\in N_{jk}}{\Gamma }_j^{k,i}}.$$

$${\mathcal{S}}_{j,k}:=\overline{\bigcup _{i\in N_{jk}}{\Gamma }_k^{j,i}},{\mathcal{S}}_{j,k}^{\mathrm{int}}:={\mathcal{S}}_{j,k}\setminus \partial {\Omega }_j,$$

$$\bigcup _{j\in N_k}{\mathcal{S}}_{j,k}={\Gamma }_k^{\mathrm{int}}.$$

$${\mathcal{S}}_j:=\bigcup _{k\in N_j}{\mathcal{S}}_{j,k},{\mathcal{S}}_j^{\mathrm{int}}:=\bigcup _{k\in N_j}{\mathcal{S}}_{j,k}^{\mathrm{int}},$$

$$\Vert u_j{\Vert }_{\infty }:=\text{ess}\underset{{\mathcal{S}}_j}{\sup }|u_j|.$$

$$\Vert u{\Vert }_{\infty }:=\underset{j=1,\dots ,N}{\max }\text{ess}\underset{{\mathcal{S}}_j}{\sup }|u_j|.$$

$$\begin{split}\mathcal{V}_{n}&:=\left\{\mathbb{v}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\operatorname{int}}_{j})\colon\|\mathbb{v}\|_{\infty}\leq 1\text{ and }\begin{cases}\mathbb{v}_{j}(x)<1\;&\forall x\in\mathcal{S}_{j}^{\operatorname{int}},\forall j\in\cup_{k=1}^{n}\mathcal{L}_{k},\\ \mathbb{v}_{j}=1&\leavevmode\nobreak\ \quad\text{on }\mathcal{S}_{j},\forall j\in\cup_{k=n+1}^{N_{\max}}\mathcal{L}_{k}.\end{cases}\right\},\\ \mathcal{C}_{n}&:=\Big{\{}\mathbb{v}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\operatorname{int}}_{j})\colon\|\mathbb{v}\|_{\infty}\leq 1\text{ and }\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{v}_{j}<1,\leavevmode\nobreak\ \forall j\in\cup_{k=1}^{n}\mathcal{L}_{k}\Big{\}}.\end{split}$$

$${\chi }_j^k:=\{\begin{array}{ll}1& \text{on }\overline{{\Gamma }_j^{k,0}}\setminus {\Gamma }_j^{\mathrm{ext}},\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\\ \in [0,1]& \text{on }\overline{{\Gamma }_j^{k,i}}\setminus {\Gamma }_j^{\mathrm{ext}}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{for }i\in N_{j,k}\text{ with }i\ne 0,\\ 0& \text{otherwise},\end{array}$$

$$\sum_{k\in{N}_{j}}\chi_{j}^{k}(x)=1\quad\text{ for all }x\in\text{int}\big{(}\Gamma_{j}^{{\operatorname{int}}}\big{)}.$$

$$\sum _{k\in N_j}{\chi }_j^k(x)=0\text{ for all }x\in {\Gamma }_j^{\text{ext}}.$$

$$\Delta w_j$$

$$w_j$$

$${\mathcal{R}}_j:u\mapsto {\mathcal{R}}_j(u)=u{|}_{{\mathcal{S}}_j},$$

$$\Delta u_j$$

$$u_j$$

$$\Delta u_j^{n+1}$$

$$u_j^{n+1}$$

$$\Delta e_j^{n+1}$$

$$e_j^{n+1}$$

$$e^n:=\left[\begin{array}{c}{e}_1^n\\ e_2^n\\ ⋮\\ e_N^n\end{array}\right].$$

$$e_j^n:=e_j^n{|}_{{\mathcal{S}}_j}.$$

$$T_{ij}:=\{\begin{array}{ll}{P}_{ij}& \text{if }j\in N_i,\\ 0& \text{otherwise. }\end{array}$$

$$P_{ij}\mathbb{v}=\mathcal{R}_{i}\Bigg{(}\mathcal{E}_{i}\Big{(}{v}_{j,i}\chi_{i}^{j}|_{\mathcal{S}_{j,i}}\Big{)}\Bigg{)}.$$

$$T=\left[\begin{array}{ccccccc}0& P_{12}& 0& 0& 0& P_{16}& P_{17}\\ P_{21}& 0& P_{23}& 0& 0& 0& P_{27}\\ 0& P_{32}& 0& P_{34}& 0& 0& P_{37}\\ 0& 0& P_{43}& 0& P_{45}& 0& P_{47}\\ 0& 0& 0& P_{54}& 0& P_{56}& P_{57}\\ P_{61}& 0& 0& 0& P_{65}& 0& P_{67}\\ P_{71}& P_{72}& P_{73}& P_{74}& P_{75}& P_{76}& 0\end{array}\right].$$

Full text

1. Introduction

Recent results in computational chemistry have raised intriguing questions about the scalability of classical one-level domain decomposition methods. In [2] the authors combined the COSMO solvation model (see, for instance [1, 15, 27]) with the classical one-level parallel Schwarz method (PSM) [17] in a sub-structured integral form. This allowed the authors to define a very efficient numerical framework capable of solving solvation problems and obtaining a fast and accurate computation of the electrostatic contribution to the solvation energy. This new framework, called ddCOSMO, exploits the physical setting of the COSMO model to define a domain decomposition. The computational domain consists of the union of van der Waals cavities associated with the atoms in the molecule, and therefore these spherical cavities can be chosen as the subdomains. The ddCOSMO implementation received a good deal of attention in the computational chemistry community (see, e.g., [20, 21]), and a major reason for its success was due to the scalability property exhibited by the classical Schwarz method.

An algorithm is said to be weakly scalable if it can solve a larger and larger problem with more and more processors in a fixed amount of time. According to classical Schwarz theory, the PSM is not scalable in general; see, e.g., [26, 3]. However, in contrast to this theory, the authors in [2] provide numerical evidence indicating that in some cases the one-level PSM converges to a given tolerance in a fixed number of iterations of a linear solver independently of the number $N$ of subdomains. This behaviour is observed if fixed-sized subdomains form a ‘chain-like’ domain and their number increases. This numerical result was subsequently rigorously proved in [4, 5, 6] for the PSM and in [3] for other one-level methods; see also [16, 13, 7] for similar scalability results.

On the other hand, it was observed by the authors that the weak scalability of the PSM is lost if the fixed-sized subdomains form a ‘globular-type’ domain $\Omega$ such that the boundaries of many subdomains lie in the interior of $\Omega$. The following question therefore arises: is it possible to quantify the lack of scalability of the PSM? To do so, for increasing $N$ one would need to estimate the number of iterations necessary to achieve a given tolerance. A typical example is to consider, for instance, the PSM for the solution of a one-dimensional Laplace problem. In [3] a heuristic argument is used to explain why for an unfortunate initialization, a contraction in the infinity norm is observed only after a number of iterations proportional to the number of subdomains $N$. A first rigorous proof of this non-scalability result was given in [8]. Interestingly, the analysis turns out to be quite complicated even in a simple one-dimensional setting.

Naturally, the problem becomes extremely challenging if one considers a three-dimensional geometry involving subdomains consisting of spherical van der Waals cavities where the simultaneous intersection of more than two subdomains is possible. To the best of our knowledge, a rigorous extension of the results given in [8] to higher dimensional problems is missing in the literature. The novelty of our work is to fill this gap and develop a general framework that allows one to study the convergence and scalability of the PSM for ‘globular-type’ domains. Our framework is based on an analysis of the infinite-dimensional Schwarz iteration operator, and as such can be viewed as a generalization of the analysis technique used in [3, 4] where finite-dimensional Schwarz iteration matrices are studied instead. Furthermore, whereas the scalability analysis in [4, 5, 6] focuses on chains of fixed sized subdomains where only simple intersections are allowed, i.e., the intersection of any three distinct subdomains is assumed to be empty, our analysis is capable of covering the case of $M$-tuple intersections, i.e., we allow an arbitrary number of subdomains to have a non-empty intersection.

It is important to remark that similar and very elegant results can be found in [12] and [22]. In [12], the authors prove that the Schwarz waveform-relaxation method applied to the heat equation contracts at most every $m+2$ iterations, $m$ being an integer representing the maximum distance of the subdomains from the boundary. However, the overlapping subdomains considered in this work are constructed by ‘artificially’ enlarging the given non-overlapping subdomains. This construction allows the authors to avoid the case that P. L. Lions refers to as “weakly overlapping” in [18] and present a convergence analysis that is based on the construction of a sequence of elliptic upper bounds, and relies on the structure of the artificially constructed overlapping domain decomposition. On the other hand, the domain decomposition considered in our paper is defined by the physics of the underlying problem and it is clearly weakly overlapping. Consequently, the proof given in [12] fails for the geometric settings considered in our work. We wish also to remark that our convergence analysis is based on a direct study of the PSM iteration operator which allows us to identify the ‘worst initialization’ and carefully track the propagation of the contraction across different subdomains comprising $\Omega$ in the course of the iterations. This is in contrast to a study of the Schwarz sequence as done in [12]. Furthermore, in order to undertake our analysis we introduce the notion of the so-called skeleton of a subdomain which allows us to carefully trace the effects of the Schwarz operator on functions defined inside the subdomains rather than on the boundaries, as done in [12].

Another interesting convergence analysis of classical Schwarz methods for several intersecting subdomains is presented in [22]. This work deals with stationary advection-reaction-diffusion problems and allows one to obtain powerful scalability results. In fact, the author’s analysis which is based on the maximum principle and so-called barrier functions even yields estimates of the contraction factor. Unfortunately, this beautiful analysis fails in the case of the Laplace equation due to the lack of advection or reaction terms, and is also not valid for domain decomposition problems in which the boundaries of many subdomains lie in the interior of the global domain. In both these cases the bounds obtained via the barrier functions become ineffective and invalidate the analysis.

The remainder of this paper is organized as follows: The main ideas and results of this work are stated in Section 2. We first introduce the notation and all the mathematical objects needed for a detailed description of the domain geometry in Subsection 2.1. Next, the PSM is formulated in Subsection 2.2, and the convergence analysis is presented in Subsection 2.3. Our main results are stated in Subsection 2.3.2, and subsequently discussed with some examples in Subsection 2.3.3. We remark that in Section 2 we restrict ourselves to a two-dimensional framework where only triple intersections are allowed. This choice is made in order to ease the notation, which is otherwise very technical and complicated, and to put more focus on the techniques used to prove our results. We show in Section 3 how to extend our analysis to the case of arbitrary types and numbers of intersections and discuss the extension to three dimensions. Section 4 contains numerical experiments that support our main results. In Section 5, we consider concrete biological molecules that have been studied in solvation models in computational chemistry (see, e.g., [19]) and understand the consequences of our analysis as it pertains to these complex bio-molecules. Finally, we present our conclusions in Section 6.

2. The Schwarz method in two dimensions

2.1. Geometric setting

We consider the Laplace equation in two dimensions. Let $\Omega\subset\mathbb{R}^{2}$ be an open, connected and bounded set, let $g\in C^{0}(\partial\Omega)$ be a given function. We must find $u\in L^{\infty}(\Omega)$ with the property that

[TABLE]

Throughout this article, we assume that the domain $\Omega$ has the following structure: Let $\Omega_{i},\leavevmode\nobreak\ i=1,\ldots,N$ be a collection of $N$ intersecting disks of radii $r_{i}>0,\leavevmode\nobreak\ i=1,\ldots,N$. Then $\Omega=\cup_{i=1}^{N}\Omega_{i}$. An example of such a geometry is shown in Figure 1. In order to avoid discussing an empty theory, we assume that $N>1$ and there exist at least two indices $i,j\in\{1,\ldots,N\}$ such that $\Omega_{i}\cap\Omega_{j}\neq\emptyset$.

We impose the following constraints on the subdomains $\Omega_{i},\leavevmode\nobreak\ i=1,\ldots,N$:

1. A1)

If two subdomains $\Omega_{i}$ and $\Omega_{j}$ have non-empty intersection, then there exist subsets of the boundaries $\Gamma_{i}\subset\partial\Omega_{i}\cap\Omega_{j}$ and $\Gamma_{j}\subset\partial\Omega_{j}\cap\Omega_{i}$ both of positive measure such that

[TABLE]

In other words, we assume that if two subdomains $\Omega_{i}$ and $\Omega_{j}$ intersect, then they must have at least a simple intersection. 2. A2)

For any choice of distinct indices $i_{1},i_{2},i_{3},i_{4}\in\{1,\ldots,N\}$ it holds that $\Omega_{i_{1}}\cap\Omega_{i_{2}}\cap\Omega_{i_{3}}\cap\Omega_{i_{4}}=\emptyset$, i.e., at most only three distinct subdomains in the collection $\{\Omega_{i}\}_{i=1}^{N}$ have non-empty intersection. In other words, we assume that our geometry consists of at most triple intersections.

Constraint A1) is imposed purely for notational convenience, and it is readily seen that the subsequent analysis does not require such an assumption. Constraint A2) is much stronger but needs to be imposed-at least initially- in order to keep the focus on the analysis rather than the notational complexities that would otherwise be introduced. In Section 3, we discuss how to weaken this constraint, and show that all of our results can be extended to the case of arbitrary types of intersections.

We now develop the necessary notation.

2.1.1. Partition of the boundary

Let $j\in\{1,\ldots,N\}$. Given the disk $\Omega_{j}$ we define the sets

[TABLE]

These sets represent a decomposition of the boundary $\partial\Omega_{j}$ into an ‘external’ part $\Gamma_{j}^{\text{ext}}$ that is common with the boundary $\partial\Omega$ of the global domain and an ‘internal’ part $\Gamma_{j}^{\operatorname{int}}$ that is contained in the global domain $\Omega$. Notice that the sets $\Gamma_{j}^{\text{ext}}$ and $\Gamma_{j}^{\text{int}}$ are both closed. An example of a domain $\Omega=\cup_{j=1}^{N}\Omega_{j}$ and the decomposition of the subdomain boundaries $\partial\Omega_{j}$ is shown in Figure 2.

Next, a further decomposition of $\Gamma_{j}^{\operatorname{int}}$ is considered. To this end, we first define the set ${N}_{j}$ as

[TABLE]

Thus, ${N}_{j}$ is the set of indices of all subdomains $\{\Omega_{i}\}_{i=1,i\neq j}^{N}$ that intersect the subdomain $\Omega_{j}$. Informally, we refer to ${N}_{j}$ as the index set of neighbours of $\Omega_{j}$. Furthermore, given the set ${N}_{j}$ we define for each $k\in{N}_{j}$ the set $N_{jk}\subset\mathbb{N}_{0}$ as

[TABLE]

Next, for each $j\in\{1,\ldots,N\}$, each $k\in{N}_{j}$ and each $i\in N_{jk}$ we define the set $\Gamma_{j}^{k,i}\subset\Gamma_{j}^{{\operatorname{int}}}$ as

[TABLE]

where $\text{int}(\cdot)$ denotes the interior of a set. Intuitively, for a fixed $k\in{N}_{j}$ and $0\neq i\in N_{jk}$, the set $\Gamma_{j}^{k,i}$ denotes the portion of $\partial\Omega_{j}$ that intersects both subdomain $k$ and subdomain $i$, while $\Gamma_{j}^{k,0}$ denotes the portion of $\partial\Omega_{j}$ that intersects only subdomain $k$. We remark that the set $\Gamma_{j}^{k,i}$ is open for every $j\in\{1,\ldots,N\},k\in{N}_{j}$ and $i\in N_{jk}$.

These definitions allow us to partition the ‘interior’ part of the boundary $\Gamma_{j}^{\operatorname{int}}$ in a natural manner. Indeed we obtain that

[TABLE]

We remark that by definition, for a fixed $j\in\{1,\ldots,N\}$ and $i,k\neq 0$ it holds that $\Gamma_{j}^{k,i}=\Gamma_{j}^{i,k}$.

As the last step, we wish to introduce the notion of so-called skeletons associated with each subdomain. To this end, let $j\in\{1,\ldots,N\}$ and $k\in{N}_{j}$ be fixed. Then we define the sets $\mathcal{S}_{j,k}^{{\operatorname{int}}}$ and $\mathcal{S}_{j,k}$ as

[TABLE]

Intuitively, for a fixed subdomain $\Omega_{j}$ and its neighbouring subdomain $\Omega_{k}$, the set $\mathcal{S}_{j,k}^{\operatorname{int}}$ consists of the closure of all the interior boundaries of the subdomain $\Omega_{k}$ that intersect $\Omega_{j}$ excluding points on the boundary $\partial\Omega_{j}$. Similarly, the set $\mathcal{S}_{j,k}$ consists of the closure of all the interior boundaries of the subdomain $\Omega_{k}$ that intersect $\overline{\Omega}_{j}$. It is now easy to see that for a fixed subdomain $\Omega_{j}$ it holds that

[TABLE]

Finally, for a fixed $j\in\{1,\ldots,N\}$ we define the sets $\mathcal{S}_{j}$ and $\mathcal{S}_{j}^{{\operatorname{int}}}$ as

[TABLE]

and we say that $\mathcal{S}_{j}$ and $\mathcal{S}_{j}^{{\operatorname{int}}}$ are the skeleton and the interior skeleton of the subdomain $\Omega_{j}$. An example of the skeleton and interior skeleton of a subdomain $\Omega_{j}$ is given in Figure 3.

Notation: Let u_{j}\in L^{\infty}\big{(}\mathcal{S}_{j}\big{)} be a function. Then we define the infinity norm of $u_{j}$, denoted $\|u_{j}\|_{\infty}$, as

[TABLE]

Additionally, let {u}_{j}\in L^{\infty}\big{(}\mathcal{S}_{j}\big{)} for each $j=1,\ldots,N$ be a family of functions and define $\mathbb{u}:=(u_{1},u_{2},\ldots,u_{N})$. Then we define the infinity norm of $\mathbb{u}$, denoted $\|\mathbb{u}\|_{\infty}$ as

[TABLE]

Finally, we remark that if $\mathbb{u}=(u_{1},u_{2},\ldots,u_{N})$ is a continuous function, i.e., each $u_{j},\leavevmode\nobreak\ j=1,\ldots,N$ is continuous on $\mathcal{S}_{j}$ then we write $\mathbb{u}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}_{j})$. Similarly, if each $u_{j},\leavevmode\nobreak\ j=1,\ldots,N$ is continuous on $\mathcal{S}^{\text{int}}_{j}$ then we write $\mathbb{u}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\text{int}}_{j})$.

2.1.2. Graph and layers of a domain

Definition 2.1** (Graph of a domain $\Omega=\cup_{j=1}^{N}\Omega_{j}$).**

Consider a domain $\Omega=\cup_{j=1}^{N}\Omega_{j}$. We define the undirected graph $\mathcal{G}$ associated with the domain $\Omega$ as the set of vertices $\mathcal{V}$ given by $\mathcal{V}:=\{1,\ldots,N\}$, and the set of edges $E$ given by $E:=\left\{(i,j)\in\mathcal{V}\times\mathcal{V}\colon j\in N_{i}\right\}$. Furthermore, we say that the node $i\in\mathcal{V}$ is a boundary node if $\partial\Omega_{i}\cap\partial\Omega$ is a set of measure greater than zero.

The graph $\mathcal{G}$ associated with a domain $\Omega$ provides an easy and intuitive visualization of the connectivity of the domain. As we shall see, the convergence properties of the Schwarz method depend on the connectivity of the domain $\Omega$.

Definition 2.2** (Layers of a graph).**

Given a domain $\Omega=\cup_{j=1}^{N}\Omega_{j}$, consider the graph $\mathcal{G}$ associated with this domain. Then we define the layers of the graph $\mathcal{G}$ in an inductive manner as follows:

1. (1)

Layer $1$ is the set of all boundary nodes of $\mathcal{G}$. This set is denoted by $\mathcal{L}_{1}$. 2. (2)

For any $j>1$ we define the graph $\mathcal{G}_{j}$ iteratively as the set $\mathcal{W}_{j}:=\mathcal{V}-\cup_{k=1}^{j-1}\mathcal{L}_{k}$ together with the associated set of edges. If $\mathcal{W}_{j}$ is non-empty, then Layer $j$ is the set $\mathcal{L}_{j}$ of all boundary nodes of $\mathcal{G}_{j}$.

Moreover, $N_{\max}$ denotes the total number of layers in the graph, and we say that the domain $\Omega$ has $N_{\max}$ layers.

An example of a domain and its decomposition into layers is shown in Figure 4.

Definition 2.3**.**

Let $\Omega=\cup_{j=1}^{N}\Omega_{j}$ be a domain with $N_{\max}$ layers. Then for any natural number $n\leq N_{\max}$ we define the sets $\mathcal{V}_{n}$ and $\mathcal{C}_{n}$ as follows:

[TABLE]

Intuitively, a non-negative function $\mathbb{v}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}_{j})$ is in the set $\mathcal{V}_{n}$ if $\mathbb{v}_{j}(x)<1$ for all points $x$ on the interior skeleton of all subdomains $\Omega_{j}$ that belong to the first $n$ layers of the domain. Notice that there is no constraint on the behaviour of the function $\mathbb{v}_{j}$ at the endpoints of the skeleton $\mathcal{S}_{j}$. Therefore, in some sense, the function $\mathbb{v}$ has begun to experience a contraction on all subdomains in the first $n$ layers, but it cannot yet be claimed that the function $\mathbb{v}$ has infinity norm strictly smaller than one on the first $n$ layers.

On the other hand, a non-negative function $\mathbb{w}$ in $\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\rm int}_{j})$ is in the set $\mathcal{C}_{n}$ if $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{w}_{j}<1$ on all subdomains $\Omega_{j}$ that belong to the first $n$ layers of the domain. Therefore, the function $\mathbb{w}$ has already contracted on all subdomains in the first $n$ layers and thus has infinity norm less than one on these subdomains.

2.1.3. Partition of unity, extension and restriction operators

Let $j\in\{1,\ldots,N\}$, we define for each $k\in N_{j}$ a function $\chi_{j}^{k}\colon\partial\Omega_{j}\rightarrow\mathbb{R}$, continuous on \text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}, with the property that

[TABLE]

and such that

[TABLE]

We say that $\{\chi_{j}^{k}\}_{\begin{subarray}{c}k\in{N}_{j}\end{subarray}}$ are the partition of unity functions on $\Gamma_{j}^{{\operatorname{int}}}$.

Remark 1**.**

Let us consider a general domain $\Omega=\cup_{i=1}^{N}\Omega_{j}$. If for all $j=1,\ldots,N$, no set $\Gamma_{j}^{\rm ext}$ has an isolated point, then we can readily use Equation (2) to define the partition of unity functions. In the pathological case where some set $\Gamma_{j}^{\rm ext}$ has an isolated point (see Figure 5), Equation (2) does not provide a correct function definition. In this pathological situation, we can modify Definition (2) by setting at least one partition of unity function to be non-zero at this isolated point. This modification preserves the continuity requirement on the partition of unity functions.

Remark 2**.**

Consider a subdomain $\Omega_{j}$ and the partition of unity functions $\{\chi_{j}^{k}\}_{\begin{subarray}{c}k\in{N}_{j}\end{subarray}}$. It follows from the definition of these functions that

[TABLE]

Consequently, there is a subset of the boundary $\partial\Omega_{j}$ with non-zero Lebesgue measure on which all partition of unity functions are zero if and only if $\partial\Omega_{j}\cap\partial\Omega$ is a set of measure greater than zero.

Remark 3**.**

Consider the partition of unity functions defined through Equation (2). Two natural questions arise: 1) Is the assumption of continuity that we have imposed on the partition of unity functions truly necessary? 2) Does the choice of the partition of unity functions affect the iterates of the Parallel Schwarz method and the asymptotic contraction factor? It can indeed be verified that the proofs of all lemmas in Section 2.3.1 still hold if we drop the continuity assumption on the partition of unity functions. On the other hand, numerical experiments in Section 4 reveal that different choices of the partition of unity functions can lead to quantitatively slightly different Schwarz iterates.

Next, we define our harmonic extension operators. To this end, let $\Gamma_{j}\subseteq\partial\Omega_{j}$ be a non-empty open set. Then we define the mapping $\mathcal{E}_{j}\colon L^{2}(\Gamma_{j})\rightarrow C^{0}(\Omega_{j})$ as the operator with the property that given any input function $u_{j}\in L^{2}(\Gamma_{j})$, the output function $w_{j}:=\mathcal{E}_{j}(u_{j})\in C^{0}(\Omega_{j})$ is the unique solution to the Dirichlet problem

[TABLE]

Finally, we define our restriction operators. Let $C^{0,\text{pc}}(\overline{\Omega_{j}})$ denote the set of functions on $\overline{\Omega_{j}}$ that are continuous on $\Omega_{j}$ and piecewise continuous on the boundary $\partial\Omega_{j}$. Then, we define the mapping $\mathcal{R}_{j}\colon C^{0,\text{pc}}(\overline{\Omega_{j}})\rightarrow C^{0}(\mathcal{S}_{j}^{\text{int}})\cap L^{\infty}(\mathcal{S}_{j})$ as

[TABLE]

and we say that $\mathcal{R}_{j}$ is the restriction operator on $\Omega_{j}$.

Intuitively, the map $\mathcal{R}_{j}$ takes as input a continuous function ${v}_{j}$ defined on $\Omega_{j}$ with piecewise continuous values on the boundary $\partial\Omega_{j}$ and gives as output the restriction of this function on the skeleton $\mathcal{S}_{j}$. Notice that we cannot claim a priori that the output function $\mathcal{R}_{j}({v}_{j})$ is continuous on the entire skeleton $\mathcal{S}_{j}$ since it is possible that there is a jump discontinuity at the endpoints of the skeleton $\mathcal{S}_{j}$ which lie on the boundary $\partial\Omega_{j}$.

In what follows, we will often consider the composition \mathcal{R}_{j}\big{(}\mathcal{E}_{j}(\lambda)\big{)} where $\lambda\in L^{2}(\Gamma_{j})$ and $\Gamma_{j}\subset\partial\Omega_{j}$. We remark that this composition is well defined if $\lambda$ is a piecewise continuous function. A schematic representation of the mappings $\mathcal{E}_{j}$ and $\mathcal{R}_{j}$ is given in Figure 6.

2.2. The Schwarz method and the operator formulation

We wish to apply the parallel Schwarz method (PSM) in order to obtain the solution to Equation (1). To this end, we consider the natural decomposition of the domain $\Omega$ into the $N$ subdomains $\{\Omega_{i}\}_{i=1}^{N}$. Then for each $j\in\{1,\ldots,N\}$ we find the function $u_{j}\in L^{\infty}(\Omega_{j})$ such that

[TABLE]

Intuitively and according to (4), we look for harmonic functions $u_{j}$ on each subdomain $\Omega_{j}$ that satisfy appropriate boundary conditions. On the exterior boundary we impose the true data $g$, while on the interior boundary we impose Dirichlet data, which is traced from the neighbouring subdomain solutions and multiplied by our partition of unity functions.

The equivalence between the global Laplace problem (1) and the domain decomposition problem (4) will be proved later in Section 2.3.2 Theorem 2.12.

Using (4), we obtain the following implementation of the PSM: Let $u_{0}\colon\overline{\Omega}\rightarrow\mathbb{R}$ be some continuous initialization. For each $j\in\{1,\ldots,N\}$ find a sequence of functions $\{u^{n}_{j}\}_{n\in\mathbb{N}}\in L^{\infty}(\Omega_{j})$ such that $u_{j}^{0}=u_{0}|_{\Omega_{j}}$ and for each $n\in\mathbb{N}$ it holds that

[TABLE]

In order to analyse the sequence of functions $\{u^{n}_{j}\}_{n\in\mathbb{N}}$ we consider the error equation associated with (5). To this end, we define for each $j\in\{1,\ldots,N\}$ and each $n\in\mathbb{N}$ the error functions $e_{j}^{n}:=u_{j}-u_{j}^{n}$. It follows that for each $j\in\{1,\ldots,N\}$ the sequence of error functions $\{e_{j}^{n}\}_{n\in\mathbb{N}}$ satisfies the equation

[TABLE]

The key step is now to recognize that a convergence analysis of (6) is easier to perform using the notion of the skeleton that we introduced earlier. To do so, for each $n\in\mathbb{N}$ we define an $N$-dimensional vector $\mathbb{e}^{n}$ of error functions as

[TABLE]

Here, each element $\mathbb{e}^{n}_{j}$ is defined as

[TABLE]

In other words, the entry $\mathbb{e}^{n}_{j}$ of the vector $\mathbb{e}^{n}$ is the restriction of the error function $e_{j}^{n}$ on the skeleton $\mathcal{S}_{j}$ of the subdomain $\Omega_{j}$.

Next, we introduce the iteration operator corresponding to the above choice of the error vector $\mathbb{e}^{n}$. To do so, we define the $N\times N$ matrix $T$ by setting

[TABLE]

The entries $P_{ij}$ of the iteration operator $T$ are operators defined as follows: For each fixed $i\in\{1,\ldots,N\}$ and $j\in N_{i}$ the mapping $P_{ij}\colon C^{0}(\mathcal{S}_{j,i})\rightarrow C^{0}(\mathcal{S}^{\operatorname{int}}_{i})\cap L^{\infty}(\mathcal{S}_{i})$ is a linear operator such that for all ${v}_{j,i}\in C^{0}(\mathcal{S}_{j,i})$ it holds that

[TABLE]

In other words, for a given subdomain $\Omega_{i}$ and a given neighbour $\Omega_{j}$, the mapping $P_{ij}$

1. (1)

takes as input some function ${v}_{j,i}$ defined on the skeleton $\mathcal{S}_{j,i}$, i.e., the part of the interior boundary of subdomain $\Omega_{i}$ that is contained in $\Omega_{j}$, and multiplies it with the partition of unity function $\chi_{i}^{j}$, 2. (2)

extends the function ${v}\chi_{i}^{j}|_{\mathcal{S}_{j,i}}$ harmonically inside the domain $\Omega_{i}$, 3. (3)

and then yields as output the restriction of this harmonic extension on the skeleton $\mathcal{S}_{i}$.

We remark that the definition of the iteration operator $T$ implies that it is block-sparse. Indeed, the $i^{\text{th}}$ row of $T$ contains non-zero entries exactly at columns $j\in N_{i}$.

Example 1**.**

Let us consider the situation of a domain $\Omega$ consisting of seven subdomains, i.e., $\Omega={\cup_{j=1}^{7}}\Omega_{j}$ as shown in Figure 3. In this setting the iteration operator $T$ is a $7\times 7$ matrix given by

[TABLE]

Using (7) and (8) we can rewrite (6) as follows

[TABLE]

Lemma 2.4**.**

Equation (6) is equivalent to Equation (9).

Proof.

We show that (9) can equivalently be rewritten as (6). To this end, let $n\in\mathbb{N}$. Then Equation (9) can be written as

[TABLE]

Using the definition of the operators $P_{ik}$ we obtain for all $j\in\{1,\ldots,N\}$

[TABLE]

A direct inspection now shows that Equation (10) is equivalent to Equation (6). ∎

It therefore follows that in order to analyse the convergence of the sequence of error functions $\{e^{n}_{j}\}_{n\in\mathbb{N}},\leavevmode\nobreak\ \leavevmode\nobreak\ j\in\{1,\ldots,N\}$ defined through Equation (6), we must study the convergence of the error vectors $\{\mathbb{e}^{n}\}_{n\in\mathbb{N}}$ defined through the matrix equation (9). More specifically, we must analyse the structure of $T$. This is the subject of the next Subsection 2.3.

We conclude this section by returning to the geometric constraint A3) that we imposed on the subdomains $\Omega_{i},\leavevmode\nobreak\ i=1,\ldots,N$ in Section 2.1.1. The following remark demonstrates the necessity of this constraint.

2.3. Convergence analysis

2.3.1. Technical lemmas

The goal of this section is to introduce some technical results that we use to prove our main theorems. We begin with the following first result.

Lemma 2.5**.**

Let $j\in\{1,\ldots,N\}$ and consider the subdomain $\Omega_{j}$, let $\widetilde{u}\in L^{\infty}(\partial\Omega_{j})$ be a non-negative, piecewise continuous function such that $\text{\rm ess}\sup_{\partial\Omega_{j}}\widetilde{u}\leq 1$. Assume that there exists at least one neighbouring index $k\in N_{j}$ of the subdomain $\Omega_{j}$ such that $\widetilde{u}<1\leavevmode\nobreak\ \text{\rm in }\mathcal{S}^{\rm{int}}_{k,j}$, and define the function $u^{\rm{new}}\in L^{\infty}(\mathcal{S}_{j})$ as $u^{\rm{new}}:=\mathcal{R}_{j}\left(\mathcal{E}_{j}(\widetilde{u})\right)$. Then it holds that

[TABLE]

An example of the application of Lemma 2.5 is given in Figure 7.

Proof.

The proof follows from the maximum principle for harmonic functions (see, e.g., [5]) combined with the definition of the partition of unity functions introduced above. ∎

Next, we prove a fundamental lemma concerning the norm of the iteration operator $T$.

Lemma 2.6**.**

Let $\mathbb{1}\in\Pi_{j=1}^{N}C^{0}\left(\mathcal{S}_{j}\right)$ denote an $N$-dimensional vector with the property that every $j^{\text{th}}$ entry of the vector is a function identically equal to 1 on the skeleton $\mathcal{S}_{j}$. Then for any natural number $n\in\mathbb{N}$ it holds that

[TABLE]

Proof.

We first prove an intermediate result. Let $\mathbb{v},\mathbb{u}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\rm int}_{j})$ be two functions such that

[TABLE]

We prove that

[TABLE]

Define ${\mathbb{w}}:=T(\mathbb{u}-\mathbb{v})$. Then by definition of $T$ we obtain that for each $j\in\{1,\ldots,N\}$ it holds that

[TABLE]

In other words, the function ${\mathbb{w}}_{j}$ is simply the restriction onto the skeleton $\mathcal{S}_{j}$ of the harmonic solution to a Dirichlet problem with boundary data from $\mathbb{u}-\mathbb{v}$.

Recall that we have by hypothesis that for each $j\in\{1,\ldots,N\}$ it holds that $\mathbb{u}_{j}-\mathbb{v}_{j}\geq 0$ on $\mathcal{S}_{j}$. It follows from the maximum principle that for each $j\in\{1,\ldots,N\}$ we have

[TABLE]

Therefore, we obtain that

[TABLE]

This completes the proof for Equation (11).

We now proceed to the proof of Lemma 2.6. Note that it suffices to prove that for all natural numbers $n\in\mathbb{N}$ and all functions $\mathbb{v}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\rm int}_{j})$ such that $\|\mathbb{v}\|_{\infty}=1$ it holds that

[TABLE]

The proof follows easily by induction from Equation (11). Indeed, the base case follows immediately by picking $\mathbb{u}=\mathbb{1}$. Next, assume that (12) holds for some $k\in\mathbb{N}$. We define functions $\widetilde{\mathbb{v}}:=T^{k}\mathbb{v}$ and $\widetilde{\mathbb{u}}:=T^{k}\mathbb{1}$ and recognize that the induction hypothesis implies that

[TABLE]

Therefore, applying once again Equation (11) to the functions $\widetilde{\mathbb{v}}$ and $\widetilde{\mathbb{u}}$ yields the required result. This completes the proof. ∎

We are now ready to state our second fundamental result.

Lemma 2.7**.**

Let $\Omega=\cup_{j=1}^{N}\Omega_{j}$ be a domain with $N_{\max}$ layers. Then for all natural numbers $n\leq N_{\max}$ it holds that $T^{n}\mathbb{1}\in\mathcal{V}_{n}$.

Proof.

The proof proceeds by induction. We consider the base case $n=1$ and show that $T\mathbb{1}\in\mathcal{V}_{1}$. To this end, let $\mathbb{v}=T\mathbb{1}$. We first consider a subdomain $\Omega_{j}$ such that the node $j$ is in $\mathcal{L}_{1}$. Clearly, we have that $\partial\Omega_{j}\cap\partial\Omega\neq\emptyset$. Therefore applying Lemma 2.5 yields that $\mathbb{v}_{j}(x)<1$ for all $x\in\mathcal{S}_{j}^{\operatorname{int}}$. Next, consider a subdomain $\Omega_{j}$ such that node $j$ is in $\mathcal{L}_{2}$ or higher. It follows from the definition of the layers that $\Omega_{j}$ does not correspond to a boundary node, and thus $\partial\Omega_{j}\cap\partial\Omega=\emptyset$. Using the definition of the iteration operator $T$, Equation (3), and Remark 2 we have that

[TABLE]

We therefore obtain that ${\mathbb{v}}_{j}=(\mathbb{1})_{j}$. This completes the proof for the base case.

Assume now that the result holds for some natural number $n=k<N_{\max}$, i.e., that $T^{k}\mathbb{1}\in\mathcal{V}_{k}$. We must prove that $T^{k+1}\mathbb{1}\in\mathcal{V}_{k+1}$.

Let $\mathbb{v}:=T^{k}\mathbb{1}$, let $\mathbb{w}:=T^{k+1}\mathbb{1}$ and consider an arbitrary subdomain $\Omega_{j}$. If node $j$ is in $\mathcal{L}_{1}$, then Lemma 2.5 yields that $\mathbb{v}_{j}(x)<1$ for all $x\in\mathcal{S}_{j}^{\operatorname{int}}$. On the other hand, if node $j$ is in $\mathcal{L}_{m}$ where $m\in\{2,\ldots,k+1\}$, then it follows from the definition of the layers that there must exist a neighbouring node $i^{\prime}\in N_{j}$ such that node $i^{\prime}$ is in $\mathcal{L}_{m-1}$. The definition of the iteration operator $T$ therefore yields

[TABLE]

The induction hypothesis implies that $\mathbb{v}_{i^{\prime}}(x)<1$ for all $x\in\mathcal{S}_{i^{\prime}}^{\operatorname{int}}$. This implies in particular that $\mathbb{v}_{i^{\prime}}(x)<1$ for all $x\in\mathcal{S}_{j,i^{\prime}}^{\operatorname{int}}$. Since the function $\mathbb{w}_{j}$ is simply the restriction on to the skeleton $\mathcal{S}_{j}$ of the harmonic extension of boundary data from the function $\mathbb{v}$, we can apply Lemma 2.5 to obtain that $\mathbb{w}_{j}(x)<1$ for all $x\in\mathcal{S}_{j}^{{\operatorname{int}}}$.

If $k+1=N_{\max}$ then we are done. If not, then consider a subdomain $\Omega_{j}$ such that the node $j$ is in $\mathcal{L}_{\widetilde{m}}$ where $\widetilde{m}\in\{k+2,\ldots,N_{\max}\}$. We must show that $\mathbb{w}_{j}=1$ on $\mathcal{S}_{j}$.

It follows from the definition of the layers that all neighbouring nodes $\ell\in N_{j}$ must belong to $\mathcal{L}_{m^{\prime}}$ where $m^{\prime}\in\{k+1,\ldots,N_{\max}\}$. Therefore, the induction hypothesis implies that for each $\ell\in N_{j}$ it holds that $\mathbb{v}_{\ell}=1$ on $\mathcal{S}_{j,\ell}$.

Using the definition of the iteration operator $T$ and Remark 2, and proceeding in the same manner as in Equation (13) we obtain that ${\mathbb{w}_{j}}=(\mathbb{1})_{j}$, and we have thus shown that $T^{k+1}\mathbb{1}=\mathbb{w}\in\mathcal{V}_{k+1}$. This completes the proof by induction. ∎

Our next goal is to obtain an analogous result for the set $\mathcal{C}_{n}$. To this end, we require the use of the following key lemma.

Lemma 2.8**.**

Let $\Omega=\cup_{i=1}^{N}\Omega_{i}$ be a domain with $N_{\max}$ layers, let \mathbb{v}\in\Pi_{j=1}^{N}C^{0}\big{(}\mathcal{S}^{\rm int}_{j}\big{)} be such that $\|\mathbb{v}_{j}\|_{\infty}\leq 1$, and let $\mathbb{w}=T\mathbb{v}$. Consider a subdomain $\Omega_{j}$ such that for all neighbouring indices $i\in N_{j}$ it holds that $\mathbb{v}_{i}(x)<1$ for all $x\in\mathcal{S}_{i}^{\operatorname{int}}$. Then it holds that

[TABLE]

Proof.

By the definition of the iteration operator we have that

[TABLE]

Therefore, we define the function $h\colon\partial\Omega_{j}\rightarrow\mathbb{R}$ as

[TABLE]

Next, we observe that the definition of the partition of unity functions and the assumptions of Lemma 2.8 imply two key properties of the function $h$.

1. (i)

It holds that

[TABLE] 2. (ii)

It holds that $h$ is a continuous function on \text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}.

As a first step, we prove a third key property of the function $h$:

1. (iii)

It holds that $h(x)<1$ for all x\in\text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}.

To this end, let x\in\text{int}(\Gamma_{j}^{\operatorname{int}})=\text{int}\Big{(}\overline{\cup_{k\in N_{j}}\cup_{i\in N_{jk}}\Gamma_{j}^{k,i}}\Big{)}. We distinguish two cases:

1. (1)

$x\in\overline{\Gamma_{j}^{k,0}}$ for some neighbouring index $k\in N_{j}$. We again have two cases:

• •

$x\in\Gamma_{j}^{k,0}$. Recalling the definition of the interior skeleton, we obtain also that $x\in\mathcal{S}_{k,j}^{\operatorname{int}}\subset\mathcal{S}_{k}^{\operatorname{int}}$. Thus it holds that $h(x)=\mathbb{v}_{k}(x)\chi_{j}^{k}(x)$. We recall that by the definition of the partition of unity functions, it holds that $\chi_{j}^{k}(x)=1$ for $x\in\Gamma_{j}^{k,0}$. Furthermore, we have by assumption that $\mathbb{v}_{k}(x)<1$ for $x\in\mathcal{S}_{k,j}^{\operatorname{int}}\subset\mathcal{S}_{k}^{\operatorname{int}}$. We therefore conclude that $h(x)<1$.

• •

$x\notin\Gamma_{j}^{k,0}$, i.e., $x$ is a boundary point of the closed set $\overline{\Gamma_{j}^{k,0}}$. Now, either $x\in{\Gamma_{j}^{\operatorname{ext}}}$ or there exists some neighbouring index $\ell\in N_{jk}$ such that $x\in\overline{\Gamma_{j}^{k,\ell}}$. In other words, there are exactly two possibilities: either $x$ is a boundary point of the exterior boundary of $\Omega_{j}$ or $x$ is a boundary point of some triple intersection. Since $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$, we have excluded the first case. The second case $x\in\overline{\Gamma_{j}^{k,\ell}}$ is covered below. 2. (2)

$x\in\overline{\Gamma_{j}^{k,\ell}}$ for some neighbouring indices $k,\ell\in N_{j}$. Recalling once again the definition of the skeletons, we obtain that $x\in\mathcal{S}_{k,j}\subset\mathcal{S}_{k}$ and $x\in\mathcal{S}_{\ell,j}\subset\mathcal{S}_{\ell}$. It therefore follows that

[TABLE]

We again have two cases:

• •

$x\in\mathcal{S}_{k}^{\operatorname{int}}$. We have by assumption that $\mathbb{v}_{k}(x)<1$ for $x\in\mathcal{S}_{k}^{\operatorname{int}}$. On the other hand, we also know that $\mathbb{v}_{\ell}(x)\leq 1$ for $x\in\mathcal{S}_{\ell}$. Therefore, Equation (14) implies that $h(x)<1$.

• •

$x\notin\mathcal{S}_{k}^{\operatorname{int}}$. Since $x\in\mathcal{S}_{k}$, we must have that $x\in\partial\Omega_{k}$. It is readily seen that this in turn implies that $x\in\mathcal{S}_{\ell}^{\operatorname{int}}$. Therefore, we obtain that $\mathbb{v}_{\ell}(x)<1$, and $\mathbb{v}_{k}(x)\leq 1$. Hence, Equation (14) again implies that $h(x)<1$.

We conclude that $h(x)<1$ for all $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$ and therefore Property (iii) of the function $h$ also holds.

Consider now the skeleton $\mathcal{S}_{j}$ and let $\mathcal{S}_{j,k},\leavevmode\nobreak\ k\in N_{j}$ be an arbitrary arc of the skeleton. We must show that $\text{ess}\sup_{\mathcal{S}_{j,k}}\mathbb{w}_{j}<1$. This is a slightly delicate argument since the function $\mathbb{w}_{j}$ need not be continuous on the closed set $\mathcal{S}_{j,k}$ due to possible jump discontinuities at the endpoints. We therefore proceed in two steps:

1. (1)

First we show that $\mathbb{w}_{j}(x)<1$ for all $x\in\mathcal{S}_{j,k}^{\operatorname{int}}$, i.e., that the function $\mathbb{w}_{j}$ is strictly smaller than one in the interior of the skeleton arc $\mathcal{S}_{j,k}$. Property (iii) of the function $h$ yields that $h(x)<1\leavevmode\nobreak\ \forall x\in\Gamma_{j}^{\operatorname{int}}\supset\mathcal{S}_{k,j}$. Therefore Lemma 2.5 yields that $\mathbb{w}_{j}(x)<1$ for all $x\in\mathcal{S}_{j}^{\operatorname{int}}\supset\mathcal{S}_{j,k}^{\operatorname{int}}$. 2. (2)

Next we show that $\lim_{\begin{subarray}{c}x\in\mathcal{S}_{j,k}\\ x\to\partial\Omega_{j}\end{subarray}}\mathbb{w}_{j}(x)<1$. In other words we must show that the limit of the function $\mathbb{w}_{j}$ along the skeleton arc $\mathcal{S}_{j,k}$, as one approaches the endpoints is strictly smaller than one. We emphasise that this step is necessary since the function $\mathbb{w}_{j}$, a priori, may contain a jump discontinuity at the endpoints of the skeleton ***This situation can arise precisely when the subdomains weakly overlap (see [18]).

To this end, let $\hat{x}\in\partial\Omega_{j}$ denote any endpoint of the skeleton arc $\mathcal{S}_{j,k}$. Once again we have two cases:

• •

Suppose $\hat{x}\in\Gamma_{j}^{\operatorname{ext}}$. Recall that $h=0$ on $\Gamma_{j}^{\operatorname{ext}}$ from Property (i) and $h(x)<1$ for all $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$ from Property (iii). Thus, the Schwarz lemma (see, e.g., [14, Pages 632-635], [5] and [18, Section 3]) implies that

[TABLE]

where $\alpha$ is a constant that depends on the angle at which the skeleton arc $\mathcal{S}_{j,k}$ intersects the boundary $\partial\Omega_{j}$ at the point $\hat{x}$.

• •

Suppose $\hat{x}\notin\Gamma_{j}^{\operatorname{ext}}$. Then \hat{x}\in\text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}. Property (ii) of the function $h$ implies that $h$ is a continuous function on \text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}. This in turn implies that the harmonic extension $\mathcal{E}_{j}(h)$ is continuous in a neighbourhood of the point $\hat{x}$. This yields in particular that

[TABLE]

Hence, the claim holds in both cases.

It therefore follows that $\text{ess}\sup_{\mathcal{S}_{j,k}}\mathbb{w}_{j}<1$. Since the skeleton arc $\mathcal{S}_{j,k}$ was arbitrary, we obtain that

[TABLE]

∎

Consider the setting of Lemma 2.8. The careful reader will observe that the proof of Lemma 2.8 required the use of all three key properties of the boundary data $h$. In particular, we explicitly used the fact that $h$ is a continuous function on the interior boundary $\Gamma_{j}^{\text{int}}$. The continuity of the function $h$ is itself a consequence of our earlier assumption that the partition of unity functions are continuous on the interior boundary $\Gamma_{j}^{\text{int}}$. One might therefore wonder if the proof of Lemma 2.8 still holds if the continuity of the partition of unity functions is not imposed. It turns out that this continuity assumption is not necessary to prove Lemma 2.8, and we may instead use the Schwarz lemma. However, as we demonstrate in Section 4 using some numerical examples, a choice of discontinuous partition of unity functions may lead to a quantitatively slightly worse contraction of the error at each iteration. We remark that the ddCOSMO implementation uses continuous partition of unity functions in practice.

We conclude this subsection by observing that Lemma 2.8 has the following consequence.

Lemma 2.9**.**

Let $\Omega=\cup_{j=1}^{N}\Omega_{j}$ be a domain with $N_{\max}$ layers, let $n\leq N_{{\max}}-2$ be a natural number, and let $\mathbb{u}\in\mathcal{V}_{n}$. Then it holds that

[TABLE]

Proof.

Let $\mathbb{v}=T\mathbb{u}$ and let $\mathbb{w}=T^{2}\mathbb{u}$. Since $\mathbb{u}\in\mathcal{V}_{n}$ and $\|T\|_{\rm{OP},\infty}\leq 1$, we obtain that $\|\mathbb{w}\|_{\infty}\leq\|\mathbb{u}\|_{\infty}\leq 1$. Therefore, we need only show that $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{w}_{j}<1$ for all $j\in\cup_{k=1}^{n}\mathcal{L}_{k}$.

Let $j\in\cup_{k=1}^{n}\mathcal{L}_{k}$ and let $\ell\in N_{j}$ be any neighbouring index. By the definition of the layers, we know that $\ell\in\cup_{k=1}^{n+1}\mathcal{L}_{k}$. On the other hand, since $\mathbb{v}\in\mathcal{V}_{n+1}$, we know that $\mathbb{v}_{\ell}(x)<1$ for all $x\in\mathcal{S}_{\ell}^{\operatorname{int}}$. Applying Lemma 2.8 immediately yields that $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{w}_{j}<1$. Since $j\in\cup_{k=1}^{n}\mathcal{L}_{k}$ was arbitrary, we conclude that $\mathbb{w}\in C_{n}$. ∎

2.3.2. Convergence results

We are now ready to state our main results.

Theorem 2.10**.**

Let $\Omega=\cup_{j=1}^{N}\Omega_{j}$ be a domain with $N_{\max}$ layers and let $\mathbb{u}\in\mathcal{V}_{N_{\max}}$. Then it holds that

[TABLE]

Proof.

Let $\mathbb{w}=T\mathbb{u}$. We must show that for all $j\in\{1,\ldots,N\}$ it holds that $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{w}_{j}<1$.

To this end, let $j\in\{1,\ldots,N\}$. Since $\mathbb{u}\in\mathcal{V}_{N_{\max}}$ it follows that that for all neighbouring indices $i\in N_{j}$ it holds that $\mathbb{u}_{i}(x)<1$ for all $x\in\mathcal{S}_{i}^{\operatorname{int}}$. Applying Lemma 2.8 immediately yields that $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{w}_{j}<1$. Since $j\in\{1,\ldots,N\}$ was arbitrary, this completes the proof. ∎

Theorem 2.10 has the following important consequence.

Corollary 2.11**.**

Let $\Omega=\cup_{j=1}^{N}\Omega_{j}$ be a domain with $N_{\max}$ layers. Then it holds that

[TABLE]

Proof.

Lemma 2.6 implies that $\|T^{N_{\max}+1}\|_{\rm{OP},\infty}=\|T^{N_{\max}+1}\mathbb{1}\|_{\infty}$.

Theorem 2.10 implies that $T^{N_{\max}+1}\mathbb{1}\in\mathcal{C}_{N_{\max}}$. By definition of the set $\mathcal{C}_{N_{\max}}$ we obtain that $\|T^{N_{\max}+1}\|_{\rm{OP},\infty}=\|T^{N_{\max}+1}\mathbb{1}\|_{\infty}<1$. ∎

Remark 4**.**

Consider the setting of Theorem 2.10. The relation (16) is sharp for $N_{\text{max}}>1$. See also Example 3 for the case $N_{\max}=1$.

Remark 5**.**

Consider the error equation (9). Corollary 2.11 implies that

[TABLE]

Theorem 2.10 also allows us to prove that the global Laplace problem (1) and the domain decomposition problem (4) are indeed equivalent.

Theorem 2.12**.**

Equations (1) and (4) are equivalent. Therefore, the PSM converges to the solution of the global Laplace problem (1).

Proof.

It is well known that there exists a unique solution $u\in C^{0}(\Omega)$ to Equation (1). A direct calculation then shows that the restrictions of this function $u|_{\Omega_{j}},\leavevmode\nobreak\ j\in\{1,\ldots,N\}$ on each subdomain $\Omega_{j}$ also satisfy Equation (4). Therefore, it suffices to show that Equation (4) must have a unique solution.

We argue by contradiction. For each $j\in\{1,\ldots,N\}$, let $v_{j},\widetilde{v}_{j}\in C^{0}(\Omega_{j})$ be two distinct solutions to (4) on the subdomain $\Omega_{j}$ and let $w_{j}:=v_{j}-\widetilde{v}_{j}$. It follows that

[TABLE]

Next, we define the function $\mathbb{w}\in\Pi_{j=1}^{N}C^{0}(\mathcal{S}^{\rm int}_{j})$ by taking, for each $j=1,\ldots,N$, the restriction of the function $w_{j}$ on the skeleton $\mathcal{S}_{j}$. It follows from Equation (17) and Lemma 2.4 that $\mathbb{w}$ satisfies the fixed-point equation $\mathbb{w}=T\mathbb{w}$. This implies in particular that

[TABLE]

where $N_{\max}$ is the number of layers in the domain $\Omega$. Equation (18) now yields that $\|T^{N_{\max}+1}\|_{\text{OP},\infty}\geq 1$, which contradicts Theorem 2.10. Therefore, we must have that $\mathbb{w}\equiv 0$ which implies that $v_{j}=\widetilde{v}_{j}$ for all $j=1,\ldots,N$, and thus Equation (4) has a unique solution. ∎

2.3.3. Discussion on convergence results and related examples

In this section, we discuss the convergence results obtained in Section 2.3 and use two examples to explain the heuristic behind them and demonstrate how they can be considered an extension of existing results given in the literature [5, 4]. In particular, in Example 2 we show how our results are capable of precisely tracking the propagation of the contraction across the different layers comprising the domain $\Omega$ in the course of the iterations. In Example 3, we consider a problem defined on a linear chain of collinear subdomains and demonstrate how our results can be considered an extension of existing convergence results in the literature.

Example 2**.**

We first consider a domain $\Omega$ consisting of the union of 51 subdomains grouped in 4 layers as shown in Figure 8 (left). We will describe visually the results obtained in Section 2.3 as they apply to this particular choice of domain.

To do so, we consider as initialization the function $\mathbb{u}^{0}=\mathbb{1}$ and follow the propagation of the contraction of the error through the course of the first five iterations of the PSM. We represent in grey the subdomains where the current approximation satisfies $\mathbb{u}^{n}_{j}=\mathbb{1}_{j}$. Moreover, we depict in white all the subdomains where the current approximation satisfies $\mathbb{u}^{n}_{j}(x)<1$ for all $x\in\mathcal{S}_{j}^{\operatorname{int}}$ and in red the subdomains such that $\text{ess}\sup_{\mathcal{S}_{j}}\mathbb{u}^{n}_{j}<1$. Since the initialisation is $\mathbb{u}^{0}=\mathbb{1}$, at the iteration [math] all the subdomains are grey; Figure 8 (left). After the first iteration, the current approximation is given by $\mathbb{u}^{1}=T\mathbb{u}^{0}$. Thanks to the external Dirichlet boundary condition, the subdomains in the first layer $\mathcal{L}_{1}$ have begun to experience a contraction. Hence, $\mathbb{u}^{1}\in\mathcal{V}_{1}$ and the first layer is then white, while all the other layers are still grey; see Figure 8 (right).

After the second iteration, the current approximation is $\mathbb{u}^{2}=T\mathbb{u}^{1}$. The contraction propagates and thus the subdomains in layer 2 are now white while the ones in layers 3 and 4 are still grey. Hence, $\mathbb{u}^{2}\in\mathcal{V}_{2}$ in agreement with Lemma 2.7; see Figure 9 (left).

The third iteration yields $\mathbb{u}^{3}=T\mathbb{u}^{2}$ and reveals new behaviour. On the one hand, the contraction is still propagating towards the last layer. Layer 3 has thus begun to experience a contraction and is displayed in white indicating that $\mathbb{u}^{3}\in\mathcal{V}_{3}$. On the other hand, since all the subdomains in Layer 1 were surrounded by white subdomains, Layer 1 is now displayed in red in accordance with Lemma 2.8. Hence $\mathbb{u}^{3}\in\mathcal{C}_{1}$. Only the three subdomains in layer 4 are still grey; see Figure 9 (right).

The fourth iteration yields $\mathbb{u}^{4}=T\mathbb{u}^{3}$. The contraction has propagated to Layer 4, which is displayed in white, and $\mathbb{u}^{4}\in\mathcal{V}_{4}$ in accordance with Lemma 2.7. Furthermore, all the subdomains in Layer 2 which were previously surrounded by white or red subdomains at iteration 3 are now displayed in red in accordance with Lemma 2.8. Thus, $\mathbb{u}^{4}\in\mathcal{C}_{2}$; see Figure 10 (left).

Finally, we consider the fifth iteration, which yields $\mathbb{u}^{5}=T\mathbb{u}^{4}$. Since all the subdomains in Layers 3 and 4 were white at the previous iteration, they are now displayed in red (see Figure 10 (right)) indicating that $\mathbb{u}^{5}\in\mathcal{C}_{5}$. In other words, we observe a contraction in the infinity norm after precisely $N_{\max}+1=4+1=5$ iterations. This is in agreement with Theorem 2.10.

Example 3**.**

Next, we consider a domain $\Omega$ consisting of a linear chain of collinear subdomains (unit disks) as shown in Figure 11.

The Schwarz operator $T$ corresponding to this geometry is the block-tridiagonal operator

[TABLE]

Notice how this operator resembles the Schwarz iteration operator constructed in Fourier space in the literature [3, Section 3.1]; see also [4]. In [3] the Schwarz operator is denoted by $T_{1D}$ and it has been shown that $\|T_{1D}\|_{\infty}=1$, where $\|\cdot\|_{\infty}$ is the usual infinity matrix norm. This is in agreement with Lemmas 2.6 and 2.7.

It is clear from Figure 3 that the considered chain is composed of a single layer since $\partial\Omega_{j}\cap\partial\Omega\neq\emptyset$ for all $j=1,\dots,N$. Lemma 2.7 implies that $T\mathbb{1}\in\mathcal{V}_{1}$ and Theorem 2.10 guarantees that $T^{2}\mathbb{1}\in\mathcal{C}_{1}$, which implies that one observes a contraction in the infinity norm after two iterations. This result is conservative for chains of subdomains composed of only one layer. In fact, it has been proved in [5] that the PSM for the solution of one-layer chains of subdomains contracts at each iteration; see also [6]. This result can be obtained also using our framework. Indeed, using Lemma 2.6, we can compute

[TABLE]

The definition of the operators $P_{j,k}$ implies that the functions

[TABLE]

solve the problems

[TABLE]

for each $j=1,\dots,N$. Notice that (19) is exactly given in [5, Theorem 7, formula (24)]. Following the same arguments as in [5], we notice that

[TABLE]

for any $j=1,\dots,N$, where the strict inequality follows by the maximum principle. The value $\text{ess }\sup_{\mathcal{S}_{j}}v_{j}$ can then be computed explicitly as in [5, Sections 4 and 5.1]. We therefore obtain that $\|T\|_{\text{OP},\infty}=\rho(\alpha)<1$, where $\rho(\alpha)$ is exactly the estimate of the contraction factor given in [5, Section 5.1].

We also remark that it is possible in the same way to analyse other chains of fixed-sized subdomains characterized by only one layer such as the examples given in [5, 6].

3. Extensions to more general settings

Recall that we began Section 2.1 by imposing constraints A1) and A2) on the geometry of the domain $\Omega$. The preceding analysis confirms our claim that the constraint A1) is imposed purely to assist clarity of exposition, and can be dropped by adopting minor changes in, for instance, the definition of the index set $N_{jk},j=1,\ldots,N,k\in N_{j}$.

On the other hand, the constraint A2) which limits the number of simultaneously intersecting disks to three, significantly restricts the generality of our analysis. A second important limitation of our analysis is the fact that it is, a priori, restricted to disks in two dimensions and it is not yet clear if our results can be extended to domains in three dimensions consisting of the union of intersecting balls. The goal of this section is to address these shortcomings and extend the preceding analysis.

3.1. Extension to arbitrary number of intersections

We assume the geometric setting introduced in Section 2.1 and remove the constraints A1) and A2). In other words, we have assumed that the disks $\Omega_{i},\leavevmode\nobreak\ i=1,\ldots,N$ may have any type and any number of simultaneous intersections. Our goal is now to outline step-by-step the changes that must be made to the analysis of Section 2 so that the convergence results still hold.

3.1.1. Notation

Let $j\in\{1,\ldots,N\}$. Given the disk $\Omega_{j}$ we define as before the sets

[TABLE]

The first complication arises when we attempt to define index sets characterising the number and types of intersections. In Section 2.1, our task was straightforward because we assumed that the domain $\Omega$ consisted of at most triple intersecting subdomains but we have no longer imposed this constraint. In order to deal with the current more general setting, we define these index sets in an inductive manner.

1. Step 1)

For each ${j_{0}}=1,\ldots,N$, we define the index set of neighbours ${N}_{j_{0}}$ as

[TABLE] 2. Step 2)

For each ${j_{0}}=1,\ldots,N$ and each $j_{1}\in N_{j_{0}}$, we define the index set of triple intersections as

[TABLE]

If $N_{j_{0}j_{1}}$ is an empty set for every $j_{1}\in N_{j_{0}}$ then we stop the procedure for this specific choice of $j_{0}$, and we set the natural number $M_{j_{0}}=1$. 3. Step m)

For each \big{(}j_{0},j_{1},\ldots,j_{m-1}\big{)}\in\{1,\ldots,N\}\times N_{j_{0}}\times N_{j_{0}j_{1}}\times\ldots\times N_{j_{0}j_{1}\ldots j_{m-2}} we define the index set of $m+1$-ple intersections, i.e., the set of all intersections generated by $m+1$ simultaneously intersecting subdomains as

[TABLE]

If $N_{j_{0}j_{1}\ldots j_{m-1}}$ is an empty set for every $j_{m-1}\in N_{j_{0}j_{1}\ldots j_{m-2}}$, then we stop the procedure for this specific choice of $j_{0}$, and we set the natural number $M_{j_{0}}=m-1$.

In addition, we define the constant $M=\max_{j_{0}}M_{j_{0}}$. Intuitively, the natural number $M_{j_{0}}+1$ denotes the maximum number of subdomains that simultaneously intersect the boundary of subdomain $\Omega_{j_{0}}$. Moreover, the natural number $M+1$ denotes the maximum number of simultaneously intersecting subdomains in the entire domain $\Omega$. We say that our domain consists of at most $M+1$-ple intersections, i.e., the intersections in our domain are generated by at most $M+1$ simultaneously intersecting subdomains.

The next complication arises now that we attempt to define decompositions of the interior boundary $\Gamma_{j_{0}}^{\operatorname{int}},\leavevmode\nobreak\ {j_{0}}=1,\ldots,N$. Let $j_{0}\in\{1,\ldots,N\}$. Then we define for every natural number $n\in\{1,\ldots,M_{j_{0}}\}$ the set $X^{n}_{j_{0}}$ of $M-$tuples as

[TABLE]

Intuitively, the set $X^{1}_{j_{0}}$ characterises the simple intersections of subdomain $\Omega_{j_{0}}$, the set $X^{2}_{j_{0}}$ characterises the triple intersections of subdomain $\Omega_{j_{0}}$ and so on. In particular, the set $X^{M_{j_{0}}}_{j_{0}}$ characterises the intersections of the highest order, i.e., the $M_{j_{0}}+1$-ple intersections of subdomain $\Omega_{j_{0}}$. It is natural to define also for every ${j_{0}}\in\{1,\ldots,N\}$ the set $Z_{j_{0}}:=\bigcup_{i=1}^{M_{j_{0}}}X_{j_{0}}^{i}$.

Let $j\in\{1,\ldots,N\}$. For any multi-index $\alpha\in Z_{j}$, we define the natural number $M_{\alpha}$ as

[TABLE]

In other words, $M_{\alpha}$ denotes the number of non-zero entries in a given multi-index $\alpha\in Z_{j}$. We can now define for all multi-indices $\alpha\in Z_{j}$ the set $\Gamma_{j}^{\alpha}\subset\partial\Omega_{j}$ as

[TABLE]

where we use an abuse of notation to write $\ell\notin\alpha$, which simply means that $\ell\neq\alpha_{i}$ for all $i=1,\ldots,M$.

These definitions now imply that

[TABLE]

It is now straightforward to define also the skeletons associated with the subdomain $\Omega_{j}$. Indeed, let $j\in\{1,\ldots,N\}$ and $k\in{N}_{j}$ be fixed. Then we define the sets $\mathcal{S}_{j,k}^{{\operatorname{int}}}$ and $\mathcal{S}_{j,k}$ as

[TABLE]

Similarly, we define the skeleton $\mathcal{S}_{j}$ and the interior skeleton $\mathcal{S}^{\text{int}}_{j}$ as

[TABLE]

It can be observed that with this notation

• •

The definition of the graph and layers of a domain remain unchanged.

• •

The definition of the sets $\mathcal{V}_{n}$ and $\mathcal{C}_{n}$ remains unchanged.

Of course, it becomes necessary to modify the definition of the partition of unity functions. Indeed, let $j\in\{1,\ldots,N\}$. Then we define for each $k\in N_{j}$ a function $\chi_{j}^{k}\colon\partial\Omega_{j}\rightarrow\mathbb{R}$, continuous on \text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}, with the property that

[TABLE]

and such that

[TABLE]

The harmonic extension operator $\mathcal{E}_{j}$ and the restriction operator $\mathcal{R}_{j}$ can now be defined analogously to Section 2, keeping in mind the notation we have just introduced.

3.1.2. Operator formulation and convergence analysis

Let us first consider the developments of Section 2.2. The definition of the error functions $\mathbb{e}^{n},\leavevmode\nobreak\ n\in\mathbb{N}$ in Equation 7 remains unchanged. Indeed, each element $\mathbb{e}_{j}^{n},\leavevmode\nobreak\ j=1,\ldots,N,\leavevmode\nobreak\ n\in\mathbb{N}$ is now defined as

[TABLE]

Finally, we consider the lemmas and theorems we have stated in Sections 2.3 and 2.3.2. Once again, we observe that the majority of the technical lemmas still hold under the current setting. In fact only the proof of Lemma 2.8 needs to be modified but the changes are minor and occur at only one point.

Indeed, consider the setting and proof of Lemma 2.8. We must modify the proof of Property (iii) of the boundary function $h$: i.e., the proof that $h(x)<1$ for all x\in\text{int}\big{(}\Gamma_{j}^{\operatorname{int}}\big{)}.

Let $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$. As before, we distinguish two cases.

1. (1)

Suppose $x\in\overline{\Gamma_{j}^{\alpha}}$ where $\alpha\in Z_{j}$ satisfies $\alpha_{2}=\alpha_{3},\ldots,=\alpha_{M}=0$ and $\alpha_{1}=k$ for some $k\in N_{j}$. In other words, we consider the case where $x$ belongs to the part of the boundary $\partial\Omega_{j}$ that is a simple intersection. We distinguish two cases.

• •

$x\in\Gamma_{j}^{\alpha}$. Recalling the definition of the interior skeleton, we obtain also that $x\in\mathcal{S}_{k,j}^{\operatorname{int}}\subset\mathcal{S}_{k}^{\operatorname{int}}$. Thus it holds that $h(x)=\mathbb{v}_{k}(x)\chi_{j}^{k}(x)$. Using the fact that the partition of unity functions are all bounded by one, and the fact that we have by assumption $\mathbb{v}_{k}(x)<1$ for $x\in\mathcal{S}_{k,j}^{\operatorname{int}}\subset\mathcal{S}_{k}^{\operatorname{int}}$, we conclude that $h(x)<1$.

• •

$x\notin\Gamma_{j}^{\alpha}$, i.e., $x$ is a boundary point of the closed set $\overline{\Gamma_{j}^{\alpha^{j}}}$. Now, either $x\in{\Gamma_{j}^{\operatorname{ext}}}$ or there exists some neighbouring index $\ell\in N_{jk}$ such that $x\in\overline{\Gamma_{j}^{\beta}}$ where $\beta\in Z_{j}$ satisfies $\beta_{1}=k$ and $\beta_{2}=\ell$. Since $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$, we have excluded the first case. The second case is covered below. 2. (2)

Suppose $x\in\overline{\Gamma_{j}^{\alpha}}$ where $\alpha\in Z_{j}$ satisfies $\alpha_{1}=k,\alpha_{2}=\ell$ for some $k\in N_{j}$ and some $\ell\in N_{jk}$. In other words, we consider the case where $x$ belongs to the part of the boundary $\partial\Omega_{j}$ that is a triple intersection or higher. It now holds that

[TABLE]

Notice now that $\text{for all }x\in\overline{\Gamma_{j}^{\alpha}}$ the partition of unity functions satisfy $\sum_{i=1}^{M_{\alpha}}\chi_{j}^{\alpha_{i}}(x)=1$, and we have by the assumptions of Lemma 2.8 that $\mathbb{v}_{\alpha_{i}}(x)\leq 1$ for all $i=1,\ldots,M_{\alpha}$. In fact, even more is true. Since x\in\overline{\Gamma_{j}^{\alpha^{j}}}\cap\text{int}\big{(}\Gamma_{j}^{{\operatorname{int}}}\big{)} it is not difficult to see that either $x\in\Omega_{k}$ or $x\in\Omega_{\ell}$, i.e., $x$ belongs to the interior of one of the two neighbouring subdomains $\Omega_{k}$ and $\Omega_{\ell}$. This being the case, we obtain from Lemma 2.5 that either $\mathbb{v}_{k}(x)<1$ or $\mathbb{v}_{\ell}(x)<1$. Equation (21) then implies that $h(x)<1$. This completes the claim.

We therefore conclude that $h(x)<1$ for all $x\in\text{int}(\Gamma_{j}^{\operatorname{int}})$ and therefore Property (iii) of the function $h$ holds also in the current more general setting. The remainder of the proof is identical.

It can now be readily checked that the proofs of the remaining lemma and theorems in Section 2.3 remain essentially unchanged. We conclude this sub-section by remarking that the reader can now understand why we chose to begin our analysis with certain simplifying assumptions: unfortunately, in a completely general setting, the notation turns increasingly complex and it becomes easy to lose focus on the main ideas of our analysis.

3.2. Extensions to three dimensions

We have assumed throughout this article that the domain $\Omega$ is a subset of $\mathbb{R}^{2}$ and is decomposed into subdomains $\{\Omega_{i}\}_{i=1}^{N}$ consisting of disks. It is pertinent to recall however, that the original ddCOSMO model involves solving the Laplace equation on a computational domain in $\mathbb{R}^{3}$, which is decomposed into subdomains $\{\Omega_{i}\}_{i=1}^{N}$ consisting of balls. We therefore briefly discuss the possibility of extending the preceding analysis to this three dimensional setting.

Let us first assume that the domain $\Omega\subset\mathbb{R}^{3}$ can be decomposed as the union of intersecting balls $\Omega=\cup_{i=1}^{N}\Omega_{i}$. As in Section 3.1, we assume that no two subdomains intersect at a single point. A careful study of Sections 2 and 3.1 now reveals that this three-dimensional geometric setting does not require any change in notation. In fact, the only lemma that uses the fact that each $\Omega_{i},\leavevmode\nobreak\ i=1,\ldots,N\subset\mathbb{R}^{2}$ is Lemma 2.8.

To be more concrete, assume the setting of Lemma 2.8. The final step of the proof consists of showing that

[TABLE]

In order to prove this result, we use the Schwarz lemma. This lemma is referenced by name in [14, Pages 632-635], and a detailed proof is provided in the case of two dimensions. Another proof of this lemma, which uses complex analysis and is therefore also valid only in two dimensions is given in [5]. Additionally, the Schwarz lemma is also referenced- this time in higher dimensions- by P.L. Lions in his classical paper on domain decomposition in [18][Section 3]. Lions claims that “… this follows easily from potential theory” but omits providing a proof. Unfortunately, apart from Lion’s single comment, we have been unable to find any proof for the Schwarz lemma in three dimensions. Moreover, the fact that the proof in [5] uses tools from complex analysis prevents a straightforward extension to three dimensions.

We therefore emphasise that the only bottleneck in extending the analysis of Sections 2 and 3.1 to the case of a domain $\Omega\subset\mathbb{R}^{3}$ composed of intersecting balls is the lack of a proof of the Schwarz lemma in three dimensions.

3.3. Extensions to other types of subdomains

The analysis we have presented thus far assumes that the domain $\Omega$ consists either of a union of intersecting disks in $\mathbb{R}^{2}$ or a union of intersecting balls in $\mathbb{R}^{3}$ (see Section 3.2). In this subsection, we briefly discuss the possibility of extending our analysis to domains composed of more general types of subdomains.

Indeed, the results presented in this manuscript will hold for a different choice of subdomains if

• (B1)

For this new choice of subdomains, it is possible to introduce definitions and notations analogous to the ones introduced in Section 2 such that the operator formulation of the Schwarz method as stated in Section 2.2 is well-posed.

• (B2)

For this new choice of subdomains, the technical results presented in Section 2.3 still hold.

Clearly, the well-posedness of the operator formulation of the Schwarz method as given in Section 2.2 requires that the boundary value problem (5) be well-posed on each subdomain $\Omega_{j},\leavevmode\nobreak\ j=1,\ldots,N$. This condition is satisfied if the subdomains $\{\Omega_{j}\}_{j=1}^{N}$ are, for instance, assumed to be Lipschitz. Moreover, it is not difficult to see that the definitions and notations introduced in Section 2 involving, for instance, interior and exterior boundary sets, skeletons and partition of unity functions, can easily be extended to Lipschitz subdomains $\{\Omega_{j}\}_{j=1}^{N}$. Consequently, condition (B1) holds for the case of Lipschitz subdomains.

Unfortunately, the analysis we have presented in Section 2.3 does not hold for general Lipschitz subdomains. In order to illustrate this point, let us consider the geometry displayed in Figure 12.

We define the sets $\Omega_{1}\subset\mathbb{R}^{2}$ and $\Omega_{2}\subset\mathbb{R}^{2}$ as

[TABLE]

and we define the domain $\Omega:=\Omega_{1}\cup\Omega_{2}$. Additionally, we denote by

[TABLE]

the curve that constitutes the right-side boundary of the domain $\Omega_{1}$. Let us now consider Lemma 2.8 for this choice of geometry and domain decomposition. Essentially, Lemma 2.8 claims the following: Let the function $h\colon\partial\Omega_{2}\rightarrow\mathbb{R}$ be defined as

[TABLE]

Let $w_{2}\in L^{2}(\Omega_{2})$ be the harmonic extension in $\Omega_{1}$ of $h$ and let the function $w_{1}\in L^{2}(\Omega_{1})$ be the solution to the BVP

[TABLE]

Then according to Lemma 2.8 we should have that

[TABLE]

Consider now the proof of Lemma 2.8 and let $P\in\partial\Omega_{1}$ be the point $(0,0)$, i.e., the lower endpoint of the curve $D$. The first case in Step (3) of the proof of Lemma 2.8 shows that

[TABLE]

This bound on the limit is established using the Schwarz lemma (see, e.g., [14, Pages 632-635]), which states that

[TABLE]

where $\alpha\in[0,1]$ is a constant that depends on the angle at which the curve $D$ intersects the boundary $\Gamma_{2}^{\rm int}$ at the point $P$. We now observe that the definition of the curve $D$ implies that it intersects $\Gamma_{2}^{\rm int}$ with angle zero at the point $P$ which yields $\alpha=1$. Thus, the bound (22) fails, i.e., Lemma 2.8 no longer holds. We emphasise that this situation does not arise in the case of a domain $\Omega\subset\mathbb{R}^{2}$ that is composed of the union of intersecting disks since the boundaries of disks do not intersect at zero angles.

Note that the geometry displayed in Figure 12 is known in the domain decomposition literature as an example of a domain $\Omega=\Omega_{1}\cup\Omega_{2}$ wherein the standard proof of geometric convergence of the PSM via a contraction argument fails (see, e.g., [9]). Therefore, it is not surprising that our analysis also does not hold in this case. We remark that convergence of the PSM in this case can still be proven but requires different tools (see, e.g., [18, Section 4]).

We conclude that the analysis we have presented in this manuscript does not hold for general non-convex Lipschitz subdomains.

4. Numerical experiments

The goal of this section is two-fold: First, we provide numerical evidence supporting our main convergence result Theorem 2.10, and its analogue for the case of quadruple intersecting disks. Second, we show the effect of adopting a different choice of the partition of unity functions on the convergence rates.

Figures 14(a) and 14(b) display the infinity norm of the Schwarz iterates generated through the Error Equation (9) using the initialisation $\mathbb{e}^{0}\equiv 1$ on $\Pi_{j=1}^{j=N}\mathcal{S}_{j}$ for geometries with a different number of maximum layers. The different geometries are displayed in Figure 13. We emphasise that continuous partition of unity functions were chosen for this and all subsequent computations.

Although it is slightly difficult to see this at a glance, the numerical results displayed in Figure 14(b) follow the theoretical results exactly. Indeed, the error function always displays a first contraction at iteration number $N_{\max}+1$. Furthermore, it is readily seen that the convergence rate of the error also depends on the total number of layers in the domain and degrades as $N_{\max}$ increases as shown in the log-lin plot displayed in Figure 14(a). Indeed, we observe that the slope of the plots and consequently the asymptotic contraction factor degrades as $N_{\max}$ increases.

One extremely complicated question that we did not address in the theoretical analysis of Section 2 is how to obtain a quantitive estimate of the contraction factor. Typically, one obtains an upper bound for the initial contraction factor, i.e., the initial decrease in the infinity norm of the Schwarz iterates (see, e.g., [4], [5] and [6]), which then improves considerably in the asymptotic limit. Figures 15(a) and 15(b) display the exact infinity norm of the Schwarz iterates together with the (numerically obtained) initial and asymptotic contraction factor for geometries with $N_{\max}=3$ and $N_{\max}=4$ layers. We observe immediately that the first contraction in both cases is nearly two orders of magnitude smaller than the asymptotic contraction factor. The numerics therefore suggest that even if we were able to obtain an upper bound for the initial contraction factor, that is, an estimate of $\|T^{N_{\max}+1}\|_{\rm{OP},\infty}$, it would be extremely conservative and not useful from a practical point of view.

In order to support our claim regarding the extension of Theorem 2.10 to more complicated geometries, we next repeat the above numerical computations for a computational domain consisting of quadruple intersecting subdomains. The different geometries are displayed in Figure 16.

Figures 17(b) displays the infinity norms of the Schwarz iterates generated through the Error Equation (9) using the initialisation $\mathbb{e}^{0}\equiv 1$ on $\Pi_{j=1}^{j=N}\mathcal{S}_{j}$ for quadruple intersecting geometries with a different number of maximum layers. Once again, although it might be difficult to tell at a glance, the numerical results follows the theoretical results exactly and we see a first contraction of the infinity norm at iteration number $N_{\max}+1$.

Figure 17(a) displays the infinity norms of the Schwarz iterates as the number of iterations increases for the three geometries in a log-lin plot. Once again we observe that the asymptotic contraction factor degrades as $N_{\max}$ increases. In addition, we plot in Figures 18(a) and 18(b) the exact infinity norms of the Schwarz iterates together with the (numerically obtained) initial and asymptotic contraction factor for geometries with $N_{\max}=3$ and $N_{\max}=4$ layers. As before, the initial contraction factor can be seen to be at least an order of magnitude smaller than the asymptotic convergence factor.

We conclude this section by returning to a question posed in Remark 3: Does the choice of partition of unity functions affect the iterates of the Parallel Schwarz method and the asymptotic contraction factor? We decided to compute the infinity norm of the Schwarz iterates generated through the Error Equation (9) using the initialisation $\mathbb{e}^{0}\equiv 1$ on $\Pi_{j=1}^{j=N}\mathcal{S}_{j}$ for the geometries displayed in Figure 13 using both the continuous partition of unity functions employed previously as well as discontinuous partition of unity functions as chosen in [10][Section 2, Page 4]. The discontinuous partition of unity functions are defined as follows: Given three intersecting subdomains $\Omega_{1},\Omega_{2},\Omega_{3}$

• •

We set the partition of unity functions $\chi_{1}^{2},\chi_{1}^{3}$ on $\partial\Omega_{1}$ as $\chi_{1}^{2}=1$ on $\overline{\Gamma_{1}^{2,3}}\cup\overline{\Gamma_{1}^{2,0}}$ and $\chi_{1}^{2}=0$ otherwise, and $\chi_{1}^{3}=1$ on $\overline{\Gamma_{1}^{3,0}}$ and $\chi_{1}^{3}=0$ otherwise.

• •

We set the partition of unity functions $\chi_{2}^{3},\chi_{2}^{1}$ on $\partial\Omega_{2}$ as $\chi_{2}^{3}=1$ on $\overline{\Gamma_{2}^{1,3}}\cup\overline{\Gamma_{2}^{3,0}}$ and $\chi_{2}^{3}=0$ otherwise, and $\chi_{2}^{1}=1$ on $\overline{\Gamma_{2}^{1,0}}$ and $\chi_{2}^{1}=0$ otherwise.

• •

We set the partition of unity functions $\chi_{3}^{1},\chi_{3}^{2}$ on $\partial\Omega_{3}$ as $\chi_{3}^{1}=1$ on $\overline{\Gamma_{3}^{1,2}}\cup\overline{\Gamma_{3}^{1,0}}$ and $\chi_{3}^{1}=0$ otherwise, and $\chi_{3}^{2}=1$ on $\overline{\Gamma_{3}^{2,0}}$ and $\chi_{3}^{2}=0$ otherwise.

In other words given any subdomain $\Omega_{j}$, on the portion of the boundary $\Gamma_{j}^{k,i}\subset\partial\Omega_{j}$ that intersects simultaneously with two other subdomains $\Omega_{k}$ and $\Omega_{i}$, we take boundary data from only one neighbour. Our results are displayed in Figures 19, and they suggest that

• •

While the choice of different partition of unity functions can lead to quantitatively slightly different infinity norms of the Schwarz iterates, the asymptotic contraction factor remains unchanged;

• •

The first contraction in the infinity norm is always observed at iteration number $N_{\max}+1$ regardless of the choice of the partition of unity functions.

For more information and details on the effect of the choice of partition of unity functions on the convergence of Schwarz methods, see also [11].

5. Application to bio-molecules: van der Waals and SAS cavities

The examples considered thus far in Section 4 have been ‘toy’ problems, which are interesting from a mathematical point of view but do not have a physical origin. The goal of this section is to consider actual biological molecules that have been studied in solvation models in computational chemistry (see, e.g., [19]) and to understand the consequences of the preceding analysis as it pertains to these complex bio-molecules.

As a reference, we consider the biological molecules that were analysed in the ddCOSMO numerical study presented in [19]. Following standard practice for the COSMO implicit solvation model, these molecules are represented as the union of intersecting balls with each ball corresponding to one atom in the molecule. The computational domain thus consisted of a union of open balls, and the Laplace equation with a prescribed, smooth boundary condition was solved on this computational domain. The problem setting in [19] therefore is exactly of the form of Equation (1), except of course that the problem is posed in three dimensions. It is important to note however that we must make a choice for the value of the radius of each atomic cavity. While this is a non-trivial question in general (see, e.g., the references in [24]), the most basic convention (see, e.g., [2]) is to construct the atomic cavities using 1.1 times the radii obtained from the Universal Force Field (UFF) [25]. This is the convention we adopt as a first step.

Two representative examples of these biological molecules, Carboxylase and Hiv-1-gp41 are displayed in Figure 20. It is readily seen that the Carboxylase molecule is ‘globular’ in structure and consists of thousands of atoms, and thus thousands of subdomains in the domain decomposition setting. On the other hand, the Hiv-1-gp41 molecule is more ‘chain-like’ in structure and consists of less than 600 atoms. Based on these visualisations, one could expect the Carboxylase molecule to contain a large number of layers and the Hiv-1-gp41 molecule to contain fewer layers. As our preceding analysis shows, this could result in very different computational efficiencies of the parallel Schwarz method (PSM) for the two molecules as described by Equation (5) and implemented through the ddCOSMO numerical algorithm.

Table 1 displays the wall-time of the ddCOSMO algorithm as reported in [19], and we can observe the surprising result that the PSM (as a preconditioner) seems to work extremely efficiently for all the biological molecules considered here. In particular, we observe excellent scaling of the method with respect to the number of atoms in each molecule. Furthermore, our own numerical computations, which are shown in the second row of Table 1 and displayed in Figure 21 indicate that there is minimal degradation of the asymptotic convergence rate of the ddCOSMO algorithm for these bio-molecules. In order to explain this observation, we decided to analyse carefully the structure of these protein molecules.

Table 2 contains information about geometrical aspects of each of these protein molecules that affect the convergence rate of the PSM. We indicate the total number of atoms, the number of layers, which is the central component of our analysis, the average number of neighbours, the average over different subdomains of the maximum degree of intersection, and the average over different subdomains of the overlap distance between neighbouring atoms. The following two observations are now key:

• •

There are exactly two layers in each bio-molecule;

• •

The other geometrical parameters which could potentially affect the convergence rate of the PSM do not vary wildly across the different molecules.

These observations help explain our numerical results and also explain why the ddCOSMO numerical algorithm works so well in practice even for apparently globular biological molecules [20, 21, 2].

In the computational chemistry community, it is generally argued that using scaled van der Waals radii to construct atomic cavities leads to a sub-optimal definition of the molecular cavity that is unable to accurately capture the physical behaviour of the molecule in a polarisable medium. A more elaborate definition of the molecular cavity is given by the so-called Surface Accessible Surface (SAS) (see, e.g., [24]). According to this definition, each atomic cavity is constructed using the sum of two radii: the atomic van der Waals radius and a so-called ‘probe radius’, which can generally be taken to be 1.4 Angstrom, i.e., the approximate radius of a water molecule. As a second step therefore, we adopt the SAS convention to construct the atomic cavities. Figure 22 displays the SAS cavity-based atomic structure of the biological molecules Carboxylase and Hiv-1-gp41 (c.f., Figure 20).

We consider once again the six biological molecules analysed in [19]. Table 3 contains information about various geometrical aspects of each of these protein molecules using the SAS cavity definition. We observe immediately that the number of layers in each molecule is now greater than two while the globular molecules UBCH5B and Carboxylase have six and eight layers respectively. Furthermore, the average over different subdomains of the maximum degree of intersection and the average over different subdomains of the overlap between neighbouring atoms are both higher than before but do not differ significantly between the molecules. Thus, the only significant differences between these molecules are now the number of layers and the average number of neighbours. Taken together, these observations suggest that the ddCOSMO algorithm should display slower convergence behaviour for all molecules and worse behaviour for the globular molecules.

As a test case, we consider the more ‘chain-like’ molecules Vancomycin and Hiv-1-gp41 and the ‘globular’ molecule UBCH5B. Notice that both the Vancomycin and Hiv-1-gp41 molecules consist of 3 layers when using the SAS definition of the cavity. On the other hand, the molecule UBCH5B, due to it’s globular structure, consists of 6 layers. We solve the ddCOSMO system of equations using the PSM. Table 4 lists the number of iterations required to solve the linear system up to a fixed tolerance.

As expected, using the SAS definition of cavities results in a larger number of iterations being required to solve the underlying linear system. It is important to note however, that even though the Hiv-1-gp41 molecule contains about $40\%$ more atoms than Vancomycin, the associated linear system for Hiv-1-gp41 can be solved in a fewer number of iterations than the linear system for Vancomycin. On the other hand, the linear system associated with the UBCH5B molecule, which contains 6 layers, now requires a significantly larger number of iterations to solve. These observations suggest the following conclusions:

• •

The key geometrical quantities which determine fast or slow convergence of the ddCOSMO algorithm for a given molecule are the number of layers (not the number of atoms by themselves) and the number of neighbouring atoms.

• •

For chain-like molecules, it is possible to use the SAS definition of the cavity without affecting the number of layers too significantly, and hence the convergence properties of ddCOSMO too adversely. On the other hand, for globular molecules, using the SAS definition results in a larger number of layers and consequently worse convergence properties.

We therefore observe that even though the SAS definition might represent a more physical representation of molecular cavities, the use of this definition for globular molecules results in a very serious deterioration in the computational efficiency of ddCOSMO.

6. Conclusions

In this work, we presented a detailed convergence and scalability analysis of the PSM for the solution of the Laplace equation defined on a domain obtained as the union of several intersecting subdomains. The geometric framework we considered is capable of describing solvation problems of interests in practical applications. Our analysis, characterizing the properties of the infinite-dimensional Schwarz operator, allowed us to prove that a first contraction in the infinity norm can be achieved in at most $N_{\max}+1$ iterations. Here $N_{\max}$ is the number of layers and represents the maximal distance of the subdomains from the exterior boundary. Numerical experiments and an application to real biological molecules demonstrate the effectiveness of our analysis.