Non-Overlapping Domain Decomposition via BURA Preconditioning of the Schur Complement

Abstract

A new class of high-performance preconditioned iterative solution methods for large-scale finite element method (FEM) elliptic systems is proposed and analyzed. The non-overlapping domain decomposition (DD) naturally introduces coupling operator at the interface $\gamma$. In general, $\gamma$ is a manifold of lower dimensions. At the operator level, a key property is that the energy norm associated with the Steklov-Poincaré operator is spectrally equivalent to the Sobolev norm of index 1/2. We define the new multiplicative non-overlapping DD preconditioner by approximating the Schur complement using the best uniform rational approximation (BURA) of ${\mathbb{L}}_{\gamma }^{1/2}$. Here, ${\mathbb{L}}_{\gamma }^{1/2}$ denotes the discrete Laplacian over the interface $\gamma$. The goal of the paper is to develop a unified framework for analysis of the new class of preconditioned iterative methods. As a final result, we prove that the BURA-based non-overlapping DD preconditioner has optimal computational complexity $O\left(n\right)$, where n is the number of unknowns (degrees of freedom) of the FEM linear system. All theoretical estimates are robust, with respect to the geometry of the interface $\gamma$. Results of systematic numerical experiments are given at the end to illustrate the convergence properties of the new method, as well as the choice of the involved parameters.

1. Introduction

This paper is aimed at constructing an analysis of methods for efficiently solving the systems of linear algebraic equations resulting from FEM discretization of second-order elliptic boundary value problems on general bounded domains $\mathsf{\Omega }\subset {\mathbb{R}}^d$. Although the approach is applicable to a wider class of multidimensional problems, we restricted our presentation to the case $\mathsf{\Omega }\subset {\mathbb{R}}^2$ with a polygonal boundary $\mathsf{\Gamma }=\partial \mathsf{\Omega }$. As a model equation, we considered the Dirichlet problem

$${\displaystyle -\sum _{i,j=1}^2\frac{\partial }{\partial x_i}\left(a_{ij}\frac{\partial u}{\partial x_j}\right)=fin\mathsf{\Omega },u=0on\mathsf{\Gamma },}$$

(1)

with ${\left\{a_{ij}\right\}}_{i,j=1}^2$ uniformly positive definite, bounded, and piecewise smooth on $\mathsf{\Omega }$. The weak formulation of (1) is: Find $u\in H_0^1\left(\mathsf{\Omega }\right)$ such that

$$\mathcal{A}(u,v)=(f,v)$$

(2)

for all $v\in H_0^1\left(\mathsf{\Omega }\right):=\{v\in H^1\left(\mathsf{\Omega }\right),v{|}_{\mathsf{\Gamma }}=0\},$ where

$${\displaystyle \mathcal{A}(u,v)=\sum _{i,j=1}^2{\int }_{\mathsf{\Omega }}{a}_{ij}\frac{\partial u}{\partial x_i}\frac{\partial v}{\partial x_j}dx,(f,v)={\int }_{\mathsf{\Omega }}fvdx.}$$

Then the FEM approximation of (2) is: Find $u_h\in S_{h,0}^1\left(\mathsf{\Omega }\right)$ such that

$$\mathcal{A}(u_h,v_h)=(f,v_h)$$

(3)

for all $v_h\in S_{h,0}^1\left(\mathsf{\Omega }\right)$. Here, $S_{h,0}^1\left(\mathsf{\Omega }\right)\subset H_0^1\left(\mathsf{\Omega }\right)$ is a finite element subspace with dimension n. Thus, the numerical solution of the elliptic problem (1) is reduced to the discrete problem (3), which in turn is equivalent to the linear system.

$$\mathbb{A}\mathbf{u}=\mathbf{f}.$$

(4)

The stiffeners matrix $\mathbb{A}\in {\mathbb{R}}^{n\times n}$ is sparse, symmetric, and positive definite. Without losing the generality, we assumed that $S_{h,0}^1\left(\mathsf{\Omega }\right)$ is the subspace of piecewise linear functions (linear triangle finite elements) defined on the triangulation ${\mathcal{T}}_h\subset \mathsf{\Omega }$.

It is assumed that the system (4) is large-scale, that is $n>>1$. For such sparse linear systems, the advantages of preconditioned conjugate gradient (PCG) solution methods increase with increasing the size n. In particular, we focused on the case of non-overlapping DD preconditioning.

A key question in DD methods is how to couple the solutions in subdomains. The interface equation involves a pseudodifferential operator, commonly called the Steklov-Poincaré operator. Many of the currently used DD methods are implicit. As a general approach, they approximate the non-local interface problem by solving a suitably constructed sequence of local (at the level of subdomains) problems. To this group belongs the original Schwarz methods, as well as: Dirichlet-Neumann, Neumann-Neumann, additive and multiplicative Schwarz, Bramble-Pasciak-Schatz (BPS), finite element tearing, and interconnect (FETI) methods. For this class of problems, independence of the convergence rate (the number of PCG iterations) of discretization and decomposition parameters does not generally hold and further consideration of the problem is needed [1]. In the spirit of implicit methods, a remarkable set of earlier results has been published in a series of papers on BPS methods, starting with [2].

An alternative second approach is to address the interface problem directly. In this case, several specific issues need to be noted: (i) the Steklov-Poincaré operator is known for a very limited set of model problems with constant coefficients in simple domains; (ii) the interface problem is non-local, and the matrix obtained after its discretization is dense, which should significantly increase the computational complexity of the DD method. Some earlier results on this approach have been published in [3,4,5,6]; see also [7] and references therein. They treat simplified cases where the interface is either one line or several non-intersecting lines. Fast Fourier transform (FFT) is used there to efficiently solve the arising systems with the square root of the matrix $tridiag(-1,2,-1)$, which approximates the corresponding Schur complement.

In this paper, we propose and analyze a new non-overlapping DD method which uses the second (direct) approach mentioned above. Our results are applicable to a very general domain partitioning, including the case of interface with cross-points. The construction is based on the spectral equivalence of the energy norm associated with the related Steklov-Poincaré operator and the norm of the Sobolev space of index 1/2. More details on the discrete representation of these norms, when FEM is applied to a numerical solution of the boundary value problem, can be found in [1]. Our approach is purely algebraic and can be successfully applied to preconditioning systems with graph Laplacians. In this context, we refer to the new discrete trace theorems recently published in [8]. The first step in the considered construction is to replace (to approximate) the Schur complement $S_{\gamma }$ corresponding to the interface $\gamma$ with ${\mathbb{L}}_{\gamma }^{1/2}$, while ${\mathbb{L}}_{\gamma }$ is the discrete Laplace-Beltrami operator assembled on $\gamma$ or, in short, the discrete Laplacian over $\gamma$. According to the trace theory, the obtained preconditioner provides optimal convergence rate, regardless of the mesh parameter h of the FEM discretization. This is confirmed, for example, by the numerical tests in [1], where the spectral decomposition of ${\mathbb{L}}_{\gamma }$ has been used to solve the fractional Laplacian problems at the interface. Obviously, such a method is computationally expensive.

This brings us to the more general problem of efficient preconditioning methods for fractional Sobolev spaces. In this context, the multigrid method developed in [9] can be considered, where the abstract additive multilevel framework has been adapted. Although the fractional diffusion on the manifold $\gamma$ is nonlocal, the preconditioner has a controlled sparsity, providing optimal convergence rate if the number of coarsening steps is fixed. Another natural approach, briefly discussed in [1], is based on a truncated spectral decomposition, where a small number of vectors are used in a modified Lanczos algorithm. The presented initial results of numerical experiments feed certain expectations for further theoretical analysis of this problem-specific model basis reduction algorithm.

We propose a new non-overlapping DD preconditioner where the BURA is applied to the block ${\mathbb{L}}_{\gamma }^{1/2}$. The BURA method was originally introduced in [10]; see also the survey paper [11]. Here, we built on the results from [12], integrating the interface preconditioning technique into the DD framework. In the present paper, we use the improved version of BURA from [13].

In general, it is clear that there are different methods for spectral fractional diffusion problems which can be applied to the block ${\mathbb{A}}_{\gamma }^{1/2}$. Recently, a unified view to these methods was presented in [14], showing that all of them can be interpreted as generated by some rational approximation. Thus, BURA methods are asymptotically the best. It should also be noted that the simplified structure of the proposed algorithm is favorable for parallel implementation. In this context, it is an advantage that the rational approximation is computed offline. Among the main contributions of this paper, we noted the following: BURA preconditioners are purely algebraic, so that no regularity assumptions are made in the theoretical analysis; under standard assumptions, the proposed DD method has optimal computational complexity; the presented numerical results prove the concept of BURA based DD preconditioning.

The structure of the paper is as follows. The spectral equivalence of the non-overlapping DD preconditioners under consideration is discussed in the next section. Section 3 is devoted to the construction and basic properties of the BURA method. Then, condition number estimates of the BURA based DD preconditioners are derived in Section 4. The optimal value of the introduced weight parameter is obtained. Algorithmic issues are discussed in Section 5, proving at the end optimal computational complexity estimates. The numerical experiments presented in Section 6 have a double motivation. On the one hand, they confirm the accuracy of the derived theoretical estimates. On the other hand, they provide an additional practical insight for the efficient usage of rational approximations of a smaller degree. Finally, we give some conclusions in Section 7.

2. Non-Overlapping DD

Let the computational domain $\mathsf{\Omega }$ be partitioned into N non-overlapping subdomains ${\mathsf{\Omega }}_i$,

$$\mathsf{\Omega }=\bigcup _{i=1}^N{\mathsf{\Omega }}_i,{\mathsf{\Omega }}_i\bigcap {\mathsf{\Omega }}_j=\varnothing \forall i\ne j,$$

and let the interface $\gamma \subset \mathbb{R}$ be alighted with the FEM triangulation ${\mathcal{T}}_h\subset \mathsf{\Omega }$,

$$\gamma =\bigcup _{i=1}^N{\gamma }_i,{\gamma }_i:=\partial {\mathsf{\Omega }}_i\backslash \partial \mathsf{\Omega }.$$

We write the linear system (4) in the $2\times 2$ block form

$$\mathbb{A}\mathbf{u}=\left(\begin{array}{cc}{\mathbb{A}}_I\hfill & {\mathbb{A}}_{I\gamma }\hfill \\ {\mathbb{A}}_{\gamma I}\hfill & {\mathbb{A}}_{\gamma }\hfill \end{array}\right)\left(\begin{array}{c}{\mathbf{u}}_I\hfill \\ {\mathbf{u}}_{\gamma }\hfill \end{array}\right)=\left(\begin{array}{c}{\mathbf{f}}_I\hfill \\ {\mathbf{f}}_{\gamma }\hfill \end{array}\right)=\mathbf{f}$$

(5)

where: $n=n_I+n_{\gamma }$; ${\mathbb{A}}_I\in {\mathbb{R}}^{n_I\times n_I}$, ${\mathbf{u}}_I\in {\mathbb{R}}^{n_I}$, and ${\mathbf{f}}_I\in {\mathbb{R}}^{n_I}$ correspond to the interior mesh nodes of the subdomains; ${\mathbb{A}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }\times n_{\gamma }}$, ${\mathbf{u}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }}$, and ${\mathbf{f}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }}$ correspond to the mesh nodes over the interface $\gamma$. The blocks ${\mathbb{A}}_{\gamma I}={\left({\mathbb{A}}_{I\gamma }\right)}^T$ represent the interaction between the nodal basis functions corresponding to the interior mesh nodes of the subdomains and the interface nodes. As typical for non-overlapping DD, the block ${\mathbb{A}}_I$ has a block-diagonal form as follows:

$${\mathbb{A}}_I=\left(\begin{array}{ccccc}{\mathbb{A}}_{I,1}\hfill & & & & \\ & \ddots \hfill & & & \\ & & {\mathbb{A}}_{I,i}\hfill & & \\ & & & \ddots \hfill & \\ & & & & {\mathbb{A}}_{I,N}\hfill \end{array}\right).$$

Here, a subdomain-by-subdomain numbering of the mesh nodes (unknowns) is used. The factorization

$$\mathbb{A}=\left(\begin{array}{cc}{\mathbb{A}}_I& \\ {\mathbb{A}}_{\gamma I}& {\mathbb{S}}_{\gamma }\end{array}\right)\left(\begin{array}{cc}{\mathbb{I}}_I& {\mathbb{A}}_I^{-1}{\mathbb{A}}_{I\gamma }\\ & {\mathbb{I}}_{\gamma }\end{array}\right)$$

(6)

holds true where ${\mathbb{S}}_{\gamma }={\mathbb{A}}_{\gamma }-{\mathbb{A}}_{\gamma I}{\mathbb{A}}_{\mathbb{I}}^{-1}{\mathbb{A}}_{I\gamma }$ is the Schur complement. Now we are ready to introduce the DD preconditioner $C_{DD}$ in the form

$${\mathbb{C}}_{DD}=\left(\begin{array}{cc}{\mathbb{A}}_I& \\ {\mathbb{A}}_{\gamma I}& {\sigma }_1{\mathbb{L}}_{\gamma }^{1/2}\end{array}\right)\left(\begin{array}{cc}{\mathbb{I}}_I& {\mathbb{A}}_I^{-1}{\mathbb{A}}_{I\gamma }\\ & {\mathbb{I}}_{\gamma }\end{array}\right),$$

(7)

${\sigma }_1>0$ is a scaling parameter that will be further specified. Here, ${\mathbb{L}}_{\gamma }$ denotes the discrete Laplace-Beltrami operator assembled on $\gamma$, which is the weighted discrete Laplacian with weight equal to the trace of the coefficient matrix ${\left\{a_{ij}\right\}}_{i,j=1}^2$ on $\gamma$.

Lemma 1.

There exists positive constants $0<c_1<c_2$ that are independent of the mesh parameter h such that

$$c_1{\mathbf{v}}_{\gamma }^T{\mathbb{L}}_{\gamma }^{1/2}{\mathbf{v}}_{\gamma }\le {\mathbf{v}}_{\gamma }^T{\mathbb{S}}_{\gamma }{\mathbf{v}}_{\gamma }\le c_2{\mathbf{v}}_{\gamma }^T{\mathbb{L}}_{\gamma }^{1/2}{\mathbf{v}}_{\gamma }\forall {\mathbf{v}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }}.$$

(8)

Proof. 

The spectral equivalence estimates (8) follow, for example, from the analysis in [1], where a more general diffusion reaction problem is considered. □

Lemma 2.

The following condition number estimate holds for the non-overlapping DD preconditioner (7)

$$\kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le \frac{max\{1,{\displaystyle \frac{c_2}{{\sigma }_1}}\}}{min\{1,{\displaystyle \frac{c_1}{{\sigma }_1}}\}}.$$

(9)

Proof. 

We rewrite (7) and (8) in the form

$$\mathbb{A}={\mathbb{U}}^T\mathbb{D}\mathbb{U}\mathrm{and}{\mathbb{C}}_{DD}={\mathbb{U}}^T{\mathbb{D}}_{{\sigma }_1}\mathbb{U},$$

where

$$\mathbb{D}=\left(\begin{array}{cc}{\mathbb{A}}_I& \\ & {\mathbb{S}}_I\end{array}\right),{\mathbb{D}}_{{\sigma }_1}=\left(\begin{array}{cc}{\mathbb{A}}_I& \\ & {\sigma }_1{\mathbb{L}}_{\gamma }^{1/2}\end{array}\right),\mathrm{and}\mathbb{U}=\left(\begin{array}{cc}{\mathbb{I}}_I& {\mathbb{A}}_I^{-1}{\mathbb{A}}_{I\gamma }\\ & {\mathbb{I}}_{\gamma }\end{array}\right).$$

Then the minimal eigenvalue of the preconditioned matrix ${\mathbb{C}}_{DD}^{-1}\mathbb{A}$ is estimated as follows:

$$\begin{array}{ccccc}{\lambda }_{min}\left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\hfill & =\hfill & {\displaystyle \underset{\mathbf{v}\ne 0}{min}\frac{(\mathbb{A}\mathbf{v},\mathbf{v})}{({\mathbb{C}}_{DD}\mathbf{v},\mathbf{v})}}\hfill & =\hfill & {\displaystyle \underset{\mathbf{v}\ne 0}{min}\frac{({\mathbb{U}}^T\mathbb{D}\mathbb{U}\mathbf{v},\mathbf{v})}{({\mathbb{U}}^T{\mathbb{D}}_{{\sigma }_1}\mathbb{U}\mathbf{v},\mathbf{v})}}\hfill \\ & & & & \\ & =\hfill & {\displaystyle \underset{\mathbf{v}\ne 0}{min}\frac{(\mathbb{D}\mathbb{U}\mathbf{v},\mathbb{U}\mathbf{v})}{({\mathbb{D}}_{{\sigma }_1}\mathbb{U}\mathbf{v},\mathbb{U}\mathbf{v})}}\hfill & =\hfill & {\displaystyle \underset{\mathbf{w}\ne 0}{min}\frac{(\mathbb{D}\mathbf{w},\mathbf{w})}{({\mathbb{D}}_{{\sigma }_1}\mathbf{w},\mathbf{w})}}\hfill \\ & & & & \\ & \ge \hfill & {\displaystyle min\{1,\frac{c_1}{{\sigma }_1}\}.}\hfill & \end{array}$$

In the same way, we derive the estimate

$${\displaystyle {\lambda }_{max}\left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le max\{1,\frac{c_2}{{\sigma }_1}\}}$$

with which the proof is completed. □

The results from Lemma 2 are summarized in

Table 1. Dependence of the estimate of $\kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)$ on the scaling parameter ${\sigma }_1$.

Interval	Condition Number Estimate
${\sigma }_1<c_1<c_2$	${\displaystyle \kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le \frac{c_2}{{\sigma }_1}>\frac{c_2}{c_1}}$
$c_1\le {\sigma }_1\le c_2$	${\displaystyle \kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le \frac{c_2}{c_1}}$
$c_1<c_2<{\sigma }_1$	${\displaystyle \kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le \frac{{\sigma }_1}{c_1}>\frac{c_2}{c_1}}$

.

Corollary 1.

The recommended interval for the scaling parameter ${\sigma }_1$ of the non-overlapping DD preconditioner ${\mathbb{C}}_{DD}$ is

$$c_1\le {\sigma }_1\le c_2.$$

(10)

In this case,

$${\displaystyle \kappa \left({\mathbb{C}}_{DD}^{-1}\mathbb{A}\right)\le \frac{c_2}{c_1}.}$$

(11)

The implementation of the DD preconditioner ${\mathbb{C}}_{DD}$ involves solving linear systems with the block ${\sigma }_1{\mathbb{L}}_{\gamma }^{1/2}$. Until recently, this was a challenge when the interface $\gamma$ contained intersecting lines. The rest of the paper is devoted to developing an effective solution to this problem using the BURA method.

3. Best Uniform Rational Approximation Preconditioning

The concept of our DD preconditioner is to replace the ${\sigma }_1{\mathbb{L}}_{\gamma }^{1/2}$ block with an appropriately chosen approximation, explicitly defined by a computationally efficient approximation of its inverse. Thus, the proposed approach directly relies on the numerical methods for spectral space-fractional problems. The survey paper [11] tracks some of the recent achievements in the numerical solution of the equation ${\mathcal{B}}^{\alpha }u=f$, $0<\alpha <1$, and is published in [11], where $\mathcal{B}$ is a symmetric positive definite (SPD) operator corresponding to second order elliptic boundary value problem in a bounded domain $\mathsf{\Omega }\subset {\mathbb{R}}^d$. The methods discussed there can be classified into the following three groups: (i) based on the Dunford-Taylor integral formula or its equivalent Balakrishnan formula; (ii) based on equivalent representation of the solution using extension to: a second order elliptic problem in $\mathsf{\Omega }\times (0,\infty )\subset {\mathbb{R}}^{d+1}$ with a local operator; a pseudo-parabolic equation in the cylinder $(x,t)\in \mathsf{\Omega }\times (0,1)$; (iii) based on spectral representation of the solution and the BURA of $z^{\alpha }$ on $[0,1]$. Although different in origin, all these methods can be interpreted as some rational approximation of ${\mathcal{B}}^{-\alpha }$. Thus, the unified approach presented and analysed in [14] provides a solid basis for comparison and evaluation. In this sense, from a mathematical point of view, the BURA methods should be the best, which is supported by various numerical tests. An example of such a comparative study is available in [11].

The BURA method was originally introduced in [10]. After that, the method was consecutively improved in several papers (for more details, see the survey in [11]). Here, we will apply the variant proposed and analyzed in [13].

Let $r_{\alpha ,k}\left(z\right)$ be the BURA of $z^{\alpha }$ in $[0,1]$ in the class of rational functions of degree k belonging to the Walsh table. This means that $r_k\left(z\right)=P_k\left(z\right)/Q_k\left(z\right)$, $P_k$ and $Q_k$ are polynomials of equal degree k. By definition, $r_{\alpha ,k}\left(z\right)$ is the minimizer $r_{\alpha ,k}\left(z\right)$ for which

$$E_{\alpha ,k}:=\underset{r_k\left(z\right)\in \mathcal{R}(k,k)}{min}\underset{z\in [0,1]}{max}|z^{\alpha }-r_k\left(z\right)|.$$

The following asymptotically sharp exponential error estimate with respect to the degree k holds true [15]:

$$E_{\alpha ,k}\le 4^{\alpha +1}sin\left(\alpha \pi \right)e^{-2\pi \sqrt{\alpha k}}.$$

(12)

Following the notation entered, the BURA approximation of ${\mathbb{B}}^{-\alpha }$ is defined as

$${\mathbb{B}}^{-\alpha }\approx {\lambda }_{1,h,\mathbb{B}}^{-\alpha }{r}_{\alpha ,k}\left({\lambda }_{1,h,\mathbb{B}}{\mathbb{B}}^{-1}\right),$$

(13)

where ${\lambda }_{1,h,\mathbb{B}}$ is the smallest eigenvalue of the SPD matrix $\mathbb{B}$.

Remark 1.

The presented construction of BURA is purely algebraic, and thus the method is applicable to arbitrary SPD matrices.

Following [12], we consider the preconditioner ${\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)$ of ${\mathbb{B}}^{\alpha }$ defined by its inverse via (14) in the form

$${\left({\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)\right)}^{-1}={\lambda }_{1,h,\mathbb{B}}^{-\alpha }{r}_{\alpha ,k}\left({\lambda }_{1,h,\mathbb{B}}{\mathbb{B}}^{-1}\right).$$

(14)

Let $\kappa \left(\mathbb{B}\right)$ denote the condition number of $\mathbb{B}$. The following lemma is then proved in [12].

Lemma 3.

Assume that the degree k is large enough so that $E_{\alpha ,k}{\kappa }^{\alpha }\left(\mathbb{B}\right)<1$. Then, the following condition number estimate for the BURA preconditioner ${\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)$, $\alpha \in (0,1)$, holds true:

$$\kappa \left({\left({\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)\right)}^{-1}{\mathbb{B}}^{\alpha }\right)\le \frac{1+E_{\alpha ,k}{\kappa }^{\alpha }\left(\mathbb{B}\right)}{1-E_{\alpha ,k}{\kappa }^{\alpha }\left(\mathbb{B}\right)}.$$

(15)

Let ${\lambda }_{i,h}$ be the eigenvalues of $\mathbb{B}$ and let ${\zeta }_{i,h}={\lambda }_{i,h}/{\lambda }_{1,h}\in \left[1,\kappa \left(\mathbb{B}\right)\right]$, for all $i=1,2,\cdots ,N$. In addition to (15), following the proof of Lemma 3, as presented in [12], we obtain the following estimates of the maximal and minimal eigenvalues of ${\left({\mathbb{C}}_{\alpha ,k}^{BURA}\right)}^{-1}{\mathbb{B}}^{\alpha }$:

$${\lambda }_{max}\left({\left({\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)\right)}^{-1}{\mathbb{B}}^{\alpha }\right)\le 1+E_{\alpha ,k}{\zeta }_{i,h}^{\alpha }\left(\mathbb{B}\right)\le 1+E_{\alpha ,k}{\kappa }^{\alpha }\left(\mathbb{B}\right),$$

(16)

$${\lambda }_{min}\left({\left({\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)\right)}^{-1}{\mathbb{B}}^{\alpha }\right)\ge 1-E_{\alpha ,k}{\zeta }_{i,h}^{\alpha }\left(\mathbb{B}\right)\ge 1-E_{\alpha ,k}{\kappa }^{\alpha }\left(\mathbb{B}\right).$$

(17)

4. Condition Number Estimates of the BURA Based DD Preconditioner

The BURA based DD preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ is defined as follows:

$${\mathbb{C}}_{DD,k}^{BURA}=\left(\begin{array}{cc}{\mathbb{A}}_I& \\ {\mathbb{A}}_{\gamma I}& {\sigma }_1{\sigma }_2{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)\end{array}\right)\left(\begin{array}{cc}{\mathbb{I}}_I& {\mathbb{A}}_I^{-1}{\mathbb{A}}_{I\gamma }\\ & {\mathbb{I}}_{\gamma }\end{array}\right)$$

(18)

where ${\sigma }_1>0$ is the scaling parameter introduced in Section 2 and ${\sigma }_2>0$ will be specified later here.

Lemma 4.

Let the assumption from Lemma 3 be valid. Then there exists positive constants ${\underline{c}}_1$ and ${\overline{c}}_2$, which are independent of the mesh parameter h,

$${\underline{c}}_1<c_1<c_2<{\overline{c}}_2,$$

(19)

such that

$${\underline{c}}_1{\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\le {\mathbf{v}}_{\gamma }^T{\mathbb{S}}_{\gamma }{\mathbf{v}}_{\gamma }\le {\overline{c}}_2{\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\forall {\mathbf{v}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }}.$$

(20)

Proof. 

Let

$${\underline{c}}_1:=c_1\left(1-E_{0.5,k}{\kappa }^{0.5}\left({\mathbb{L}}_{\gamma }\right)\right)and{\overline{c}}_2:=c_2\left(1+E_{0.5,k}{\kappa }^{0.5}\left({\mathbb{L}}_{\gamma }\right)\right).$$

(21)

Thus (19) is satisfied, where the left-had inequality follows from the assumption from Lemma 3. Combining (8) and (16), for every ${\mathbf{v}}_{\gamma }\in {\mathbb{R}}^{n_{\gamma }}$, we obtain

$$\begin{array}{ccccc}\hfill {\mathbf{v}}_{\gamma }^T{\mathbb{S}}_{\gamma }{\mathbf{v}}_{\gamma }& \le & c_2{\mathbf{v}}_{\gamma }^T{\mathbb{L}}_{\gamma }^{1/2}{\mathbf{v}}_{\gamma }\hfill & \le & c_2\left(1+E_{0.5,k}{\kappa }^{0.5}\left({\mathbb{L}}_{\gamma }\right)\right){\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\hfill \\ & & & & \\ & & & :=& {\overline{c}}_2{\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\hfill \end{array}$$

and

$$\begin{array}{ccccc}\hfill {\mathbf{v}}_{\gamma }^T{\mathbb{S}}_{\gamma }{\mathbf{v}}_{\gamma }& \ge & c_1{\mathbf{v}}_{\gamma }^T{\mathbb{L}}_{\gamma }^{1/2}{\mathbf{v}}_{\gamma }\hfill & \ge & c_1\left(1-E_{0.5,k}{\kappa }^{0.5}\left({\mathbb{L}}_{\gamma }\right)\right){\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\hfill \\ & & & & \\ & & & :=& {\underline{c}}_1{\mathbf{v}}_{\gamma }^T{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right){\mathbf{v}}_{\gamma }\hfill \end{array}$$

which completes the proof. □

Lemma 5.

The following condition number estimate holds for the BURA based non-overlapping DD preconditioner (18)

$$\kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{max\{1,{\displaystyle \frac{{\overline{c}}_2}{{\sigma }_1{\sigma }_2}}\}}{min\{1,{\displaystyle \frac{{\underline{c}}_1}{{\sigma }_1{\sigma }_2}\}}}.$$

(22)

Proof. 

The proof follows from (20), applying analogous arguments as in the proof of Lemma 2. □

The relative condition number of the BURA based DD preconditioner is controlled by the product of the positive scaling parameters ${\sigma }_1$ and ${\sigma }_2$. It is also important that the positive constants ${\underline{c}}_1$ and ${\overline{c}}_2$ are explicitly defined in (19). Lemma 5 provides an analysis of how the estimate of $\kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)$ depends on the scaling parameters. The results are summarized in

Table 2. Dependence of the estimate of $\kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)$ on the product of scaling parameters ${\sigma }_1{\sigma }_2$.

Interval	Condition Number Estimate
${\sigma }_1{\sigma }_2<{\underline{c}}_1<{\overline{c}}_2$	${\displaystyle \kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{{\overline{c}}_2}{{\sigma }_1{\sigma }_2}>\frac{{\overline{c}}_2}{{\underline{c}}_1}}$
${\underline{c}}_1\le {\sigma }_1{\sigma }_2\le {\overline{c}}_2$	${\displaystyle \kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{{\overline{c}}_2}{{\underline{c}}_1}}$
${\underline{c}}_1<{\overline{c}}_2<{\sigma }_1{\sigma }_2$	${\displaystyle \kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{{\sigma }_1{\sigma }_2}{{\underline{c}}_1}>\frac{{\overline{c}}_2}{{\underline{c}}_1}}$

.

Corollary 2.

The recommended interval for the product of the scaling parameters ${\sigma }_1{\sigma }_2$ involved in the the definition of the BURA based non-overlapping DD preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ is

$${\underline{c}}_1\le {\sigma }_1{\sigma }_2\le {\overline{c}}_2.$$

(23)

In this case,

$${\displaystyle \kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{{\overline{c}}_2}{{\underline{c}}_1}.}$$

(24)

The behaviour of the preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ for large degree k is determined by the following asymptotic properties, which follow from (12) and (23):

$$\underset{k\to \infty }{lim}{\underline{c}}_1=c_1,\underset{k\to \infty }{lim}{\overline{c}}_2=c_2,\underset{k\to \infty }{lim}{\sigma }_2=1$$

and consequently

$${\displaystyle \kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le \frac{c_2}{c_1}+{\epsilon }_k,where\underset{k\to \infty }{lim}{\epsilon }_k=0.}$$

(25)

The error of the BURA has an exponential convergence rate with respect to the degree k; see (12). However, the use of smaller k may also be of practical interest when the condition $E_{\alpha ,k}^{BURA}\left({\mathbb{L}}_h\right)<1$ is not met. Such cases are numerically studied in [12], showing that for coupled interface problems and $\alpha =0.5$, $k\ge 4$ may be enough. This conclusion is also supported by the results of the numerical experiments presented in Section 6.

5. Implementation Issues and Computational Complexity Estimate

An additive representation of the rational function ${\tilde{r}}_{\alpha ,k}\left(z\right):=r_{\alpha ,k}(1/z)$ can be used for implementation of the BURA method:

$${\tilde{r}}_{\alpha ,k}\left(z\right)={\tilde{c}}_0+\sum _{i=1}^k\frac{{\tilde{c}}_i}{z-{\tilde{d}}_i}$$

where ${\tilde{c}}_i>0$ and ${\tilde{d}}_i<0$ [13]. In this way, the inverse of BURA approximation (14) reads as

$${\left({\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)\right)}^{-1}={\lambda }_{1,h}^{-\alpha }\left({\tilde{c}}_0\mathbb{I}+\sum _{i=1}^k\left({\lambda }_{1,h}{\tilde{c}}_i\right){(\mathbb{B}-{\lambda }_{1,h}{\tilde{d}}_i\mathbb{I})}^{-1}\right).$$

(26)

Therefore, the implementation of the preconditioner ${\mathbb{C}}_{\alpha ,k}^{BURA}\left(\mathbb{B}\right)$ requires solving k auxiliary linear systems with sparse SPD matrices, which are obtained by positive diagonal perturbations of $\mathbb{B}$.

The following assumptions will be used to analyze the computational complexity of the introduced BURA based non-overlapping DD preconditioner:

(A1)

A quasi uniform mesh with a mesh parameter h is used for FEM approximation of the considered elliptic problem. Thus, the number of unknowns is $n=O\left(h^{-2}\right)$ and the number of unknowns related to the interface $\gamma$ is $n_{\gamma }=O\left(h^{-1}\right)\le c_{A1}\sqrt{n}$.

(A2)

Under the assumption A1, for the condition number of the stiffness matrix holds, $\kappa \left(\mathbb{A}\right)=O(h^{-2}$). We will assume in particular that $\kappa \left({\mathbb{L}}_{\gamma }\right)\le c_{A2}n$.

(A3)

The stiffness matrix $\mathbb{A}$ has at most $c_{A3}$ non=zero elements per row (column).

(A4)

For systems with sparse SPD matrices that appear in the implementation of ${\mathbb{C}}_{DD}^{BURA}$, a solver with optimal computational complexity is applied. In particular, we will assume that the cost of the solver is bounded from above by $c_{A4}$ arithmetic operations per unknown. Such commonly used solvers are based, for example, on multigrid or multilevel preconditioners.

(A5)

For the iterative solution of the system (4), the PCG method with preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ is used.

The computational complexity of one PCG iteration reads as

$${\mathcal{N}}_{\mathbb{C}}^{it}={\mathcal{N}}_{{\mathbb{C}}^{-1}}+{\mathcal{N}}_{\mathbb{A}}+10n+2$$

(27)

where ${\mathcal{N}}_{{\mathbb{C}}^{-1}}$ is the cost of solving a system with the preconditioner $\mathbb{C}$, and ${\mathcal{N}}_{\mathbb{A}}$ is the cost of a matrix vector multiplication with the stiffness matrix $\mathbb{A}$.

Thus, the next step is to estimate the computational complexity of solving a system with the preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$. Following (18), the complexity of one iteration can be written in the form

$${\mathcal{N}}_{{\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}}^{it}={\mathcal{N}}_{{\left({\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)\right)}^{-1}}+2{\mathcal{N}}_{{\mathbb{A}}_I^{-1}}+{\mathcal{N}}_{{\mathbb{A}}_{I\gamma }}+{\mathcal{N}}_{{\mathbb{A}}_{\gamma I}}+2n+{\mathcal{N}}_{\mathbb{A}}+10n+2$$

(28)

where ${\mathcal{N}}_{{\mathbb{A}}_I^{-1}}$ and ${\mathcal{N}}_{{\left({\mathbb{C}}_{0.5,k}^{BURA}\right)}^{-1}}$ are the costs for solving systems with the blocks ${\mathbb{A}}_I$ and ${\mathbb{C}}_{0.5,k}^{BURA}$, respectively, and ${\mathcal{N}}_{{\mathbb{A}}_{I\gamma }}$ and ${\mathcal{N}}_{{\mathbb{A}}_{\gamma I}}$ are the costs for matrix vector multiplications with ${\mathbb{A}}_{I\gamma }$ and ${\mathbb{A}}_{\gamma I}$. Then, (26) is applied to the last before term in (28), thus obtaining

$${\mathcal{N}}_{{\left({\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)\right)}^{-1}}=\sum _{i=1}^k{\mathcal{N}}_{{({\mathbb{L}}_{\gamma }-{\lambda }_{1,h}{\tilde{d}}_i{\mathbb{I}}_{\gamma })}^{-1}}+(3k+1)n_{\gamma }+2k$$

(29)

where ${\mathcal{N}}_{{({\mathbb{L}}_{\gamma }-{\lambda }_{1,h}{\tilde{d}}_i\mathbb{I})}^{-1}}$ is the cost of solving a system with the sparse SPD matrix ${\mathbb{L}}_{\gamma }-{\lambda }_{1,h}{\tilde{d}}_i\mathbb{I}$.

Combining (28) and (29) and using the assumptions (A1)–(A5), we get the estimate

$$\begin{array}{ccc}\hfill {\mathcal{N}}_{{\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}}^{it}& \le & \left(c_{A4}k+3(k+1)\right)n_{\gamma }+2k+c_{A4}n+c_{A3}n+2n+c_{A3}n+10n+2\hfill \\ \hfill & =& 2(c_{A3}+c_{A4}+6)n+\left(c_{A4}k+3k+1\right)n_{\gamma }+2(k+1)\hfill \\ \hfill & \le & 2(c_{A3}+c_{A4}+6)n+\left(c_{A4}k+3k+1\right)c_{A1}\sqrt{n}+2(k+1).\hfill \end{array}$$

(30)

Lemma 6.

Let k be the minimal integer such that

$$\sqrt{k}\ge \frac{ln\left(24\sqrt{c_{A2}n}\right)}{\sqrt{2}\pi }.$$

(31)

Then the following estimate is valid for the relative condition number of the non-overlapping DD preconditioner

$$\kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)\le 2\frac{c_1}{c_2}.$$

(32)

Proof. 

From assumption (A2) and (31), we get the inequality

$$\sqrt{k}\ge \frac{ln\left(24\sqrt{\kappa \left({\mathbb{L}}_{\gamma }\right)}\right)}{\sqrt{2}\pi },$$

which is equivalent to

$$e^{\sqrt{2}\pi \sqrt{k}}\ge 24\sqrt{\kappa \left({\mathbb{L}}_{\gamma }\right)}.$$

The last inequality can be written in the form

$$\sqrt{\kappa \left({\mathbb{L}}_{\gamma }\right)}8e^{-2\pi \sqrt{k/2}}\le \frac{1}{3}$$

and consequently

$$E_{0.5,k}{\left(\kappa \left({\mathbb{L}}_{\gamma }\right)\right)}^{0.5}\le \frac{1}{3}.$$

(33)

Combining the last estimate with the definitions (21) and the estimate (24) from Corollary 2, we get

$$\kappa \left({\mathbb{C}}_{DD,k}^{BURA}{)}^{-1}\mathbb{A}\right)\le \frac{{\displaystyle \frac{4}{3}{c}_1}}{{\displaystyle \frac{2}{3}{c}_2}},$$

which completes the proof. □

The presented analysis is summarized in the following theorem.

Theorem 1.

Let the assumptions (A1)–(A5) be valid and let

$$k:=⌈\frac{ln\left(24\sqrt{c_{A2}n}\right)}{\sqrt{2}\pi }⌉.$$

(34)

Then, for large scale problems (i.e., n is large enough) the PCG iterative solution of the system (4) with the BURA based non-overlapping DD preconditioner (18) has optimal computational complexity.

Proof. 

From Lemma 6 and the setting (34), we obtain that the number of PCG iterations with the preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ that are needed to achieve the prescribed accuracy is bounded by a constant that is independent of n. Thus, the proof of the theorem follows directly from (30), given that, for sufficiently large n, the first term dominates in the estimate of ${\mathcal{N}}_{{\mathbb{C}}_{DD,k}^{BURA}{)}^{-1}}^{it}$. □

Remark 2.

A larger degree of the BURA k can be used while maintaining the optimal computational complexity of the BURA based non-overlapping DD preconditioner. In this way, the estimate (33) can be improved to

$$E_{0.5,k}{\left(\kappa \left({\mathbb{L}}_{\gamma }\right)\right)}^{0.5}\le \frac{1}{m},$$

$m>3$ is a given integer, which reduces the constant in (32) from 2 to $(1+1/m)/(1-1/m)$.

6. Numerical Tests

The presented numerical tests are aimed at an experimental study of the preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$, introduced and analyzed in the previous sections. The case of Laplacian is considered, i.e., ${\left\{a_{ij}\right\}}_{i,j=1}^2=diag(1,1)$ in (1).

Example 1.

The numerical results are for a test problem in $\mathsf{\Omega }=(0,1)\times (0,1)$, $N=4$, and ${\mathsf{\Omega }}_1=(0,0.5)\times (0.5,1)$, ${\mathsf{\Omega }}_2=(0.5,1)\times (0.5,1)$, ${\mathsf{\Omega }}_3=(0,0.5)\times (0,0.5)$, and ${\mathsf{\Omega }}_4=(0.5,1)\times (0,0.5)$. The right-hand side is $f(x,y)=sin\left(\pi x\right)sin\left(\pi y\right)$, corresponding to the solution $u(x,y)=sin\left(\pi x\right)sin\left(\pi y\right)/\left(2{\pi }^2\right)$. All matrices are obtained by FEM discretization of the boundary value problem. Linear finite elements are used. A uniform rectangular mesh with a mesh parameter h is applied, in which each square cell is split into two triangles. The number of PCG iterations with the non-overlapping DD preconditioner ${\mathbb{C}}_{DD,k}^{BURA}$ are examined, where the parameters k and ${\sigma }_1{\sigma }_2$ are varied.

We start with a uniform rectangular mesh with mesh parameter $h=1/16$, $n_e=512$ triangle finite elements, and $n=255$ unknowns. Then, ${\ell }_r=7$ consecutive uniform refinement steps are applied, thus obtaining the finest mesh with $h=1/2048$, $n_e=$ 8,388,608 elements, and $n=$ 4,190,209 unknowns. The non-overlapping DD of ${\displaystyle \mathsf{\Omega }=\bigcup _{i=1}^4{\mathsf{\Omega }}_i}$, the initial mash, and the first mesh refinement are shown in Figure 1.

The first step in our analysis concerns the behaviour of the condition number $\kappa \left({\mathbb{L}}_{\gamma }\right)$, as shown in

Table 3. Number of finite elements $n_e$, number of unknowns n, number of unknowns corresponding to the interface $n_{\gamma }$, and condition number $\kappa \left({\mathbb{L}}_{\gamma }\right)$: ${\ell }_r$ levels of uniform mesh refinements.

${\mathit{\ell }}_{\mathit{r}}$	${\mathit{n}}_{\mathit{e}}$	n	${\mathit{n}}_{\mathit{\gamma }}$	$\mathit{\kappa }\left({\mathbb{L}}_{\mathit{\gamma }}\right)$
7	8,388,608	4,190,209	4093	2,248,960	*
6	2,097,152	1,046,529	2045	562,240	*
5	524,288	261,121	1101	140,560	*
4	131,072	65 025	509	35,140	*
3	32,768	16,129	253	8785
2	8192	3696	125	2179
1	2048	961	61	537
0	512	255	29	130

. The condition numbers corresponding to the first four meshes are computed. The results confirm the estimate $\kappa \left({\mathbb{L}}_{\gamma }\right)=O\left(h^{-2}\right)$. More precisely, we observe that, for smaller h, $\kappa \left({\mathbb{L}}_{\gamma }\right)$ is proportional to $h^{-2}$. The last four condition numbers are marked with $“*”$, indicating that they were obtained using the estimate.

As shown in Corollary 2, the condition number of the BURA based non-overlapping DD preconditioner tends to the optimal value of $c_2/c_1$ when $\kappa \left({\left[{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)\right]}^{-1}{\mathbb{L}}_{\gamma }^{0.5}\right)$ tends to 1. In addition to the theoretical estimates, we show in

Table 4. Computed estimates of $\kappa \left({\left[{\mathbb{C}}_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)\right]}^{-1}{\mathbb{L}}_{\gamma }^{0.5}\right)$, with respect to $3\le k\le 9$, for $\delta =\kappa \left({\mathbb{L}}_{\gamma }\right)\in \{{10}^5,{10}^6,{10}^7,{10}^8\}$.

k	$\mathit{\delta }={10}^5$	$\mathit{\delta }={10}^6$	$\mathit{\delta }={10}^7$	$\mathit{\delta }={10}^8$
3	1.46	3.23	9.98	31.50
4	1.38	1.46	3.29	10.19
5	1.08	1.43	1.46	3.78
6	1.03	1.08	1.46	1.66
7	1.02	1.06	1.18	1.46
8	1.01	1.03	1.08	1.34
9	1.00	1.01	1.03	1.08

the computed condition numbers of BURA preconditioners (see for more details [12]), with respect to $3\le k\le 9$, for $\delta =\kappa \left({\mathbb{L}}_{\gamma }\right)\in \{{10}^5,{10}^6,{10}^7,{10}^8\}$.

Combining the experimental data from Table 3 and Table 4, we conclude that $k>9$ should be enough to ensure an optimal relative condition number of the $C_{DD,k}^{BURA}$ preconditioner.

The results reported in

Table 5. Number of finite elements $n_e$, number of unknowns n, and number of PCG iterations $n_{it}$ with preconditioner $C_{DD,12}^{BURA}$: ${\ell }_r$ levels of uniform mesh refinements, ${\sigma }_1{\sigma }_2=2$, and stopping criteria $ϵ={10}^{-6}$.

${\mathit{\ell }}_{\mathit{r}}$	${\mathit{n}}_{\mathit{e}}$	n	${\mathit{n}}_{\mathit{it}}$
7	8,388,608	4,190,209	9
6	2,097,152	1,046,529	8
5	524,288	261,121	8
4	131,072	65,025	8
3	32,768	16,129	7

illustrate the optimal convergence rate of the BURA based non-overlapping DD method. The used degree of BURA is $k=12$. In this set of experiments, we have chosen the parameter ${\sigma }_1{\sigma }_2=2$. According to the theory, the number of iterations is independent of h and therefore of n.

In

Table 6. Number of PCG iterations $n_{it}$ with preconditioner $C_{DD,12}^{BURA}$ varying the scaling parameter ${\sigma }_1{\sigma }_2$: ${\ell }_r$ levels of uniform mesh refinements and stopping criteria $ϵ={10}^{-6}$.

${\ell }_{\mathit{r}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$
	${\mathit{\sigma }}_{\mathbf{1}}{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{1}$	${\mathit{\sigma }}_{\mathbf{1}}{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{2}$	${\mathit{\sigma }}_{\mathbf{1}}{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{3}$	${\mathit{\sigma }}_{\mathbf{1}}{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{4}$	${\mathit{\sigma }}_{\mathbf{1}}{\mathit{\sigma }}_{\mathbf{2}}=\mathbf{5}$
7	10	9	9	10	10
6	10	8	9	9	9
5	10	8	9	9	9
4	9	8	8	9	9
3	9	7	8	8	8

, we examine the influence of the scaling parameter by changing ${\sigma }_1{\sigma }_2\in \{1,2,3,4,5\}$ and keeping the degree of the BURA $k=12$. The presented results show a stable behaviour of the number of PCG iterations. The minimal values for $n_{it}$ are obtained for ${\sigma }_1{\sigma }_2=2$. This picture is in very good agreement with Lemma 5 and Consequence 2. It is also worth noting that, for $k=12$, the condition number $\kappa \left({\left({\mathbb{C}}_{DD,k}^{BURA}\right)}^{-1}\mathbb{A}\right)$ is very close to 1 (see the numerical experiments in [12]), resulting in ${\sigma }_2\approx 1$.

So far, the degree $k=12$ has been chosen to ensure that $E_{0.5,k}{\left(\kappa \left({\mathbb{L}}_{\gamma }\right)\right)}^{0.5}$ is sufficiently small.

Table 7. Number of PCG iterations $n_{it}$ with preconditioner $C_{DD,k}^{BURA}$ varying the degree k: ${\ell }_r$ levels of uniform mesh refinements, ${\sigma }_1{\sigma }_2=2$, and stopping criteria $ϵ={10}^{-6}$.

${\ell }_{\mathit{r}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$	${\mathit{n}}_{\mathit{it}}$
	$\mathit{k}=\mathbf{4}$	$\mathit{k}=\mathbf{8}$	$\mathit{k}=\mathbf{12}$
7	27	8	9
6	19	8	8
5	14	9	8
4	9	8	8
3	8	7	7

presents a comparative analysis of the number of iterations when the degree of BURA, $k\in \{4,8,12\}$, varies. We notice that $k=4$ is applicable for smaller problems, with up to $n=65,025$ unknowns, while in general $k\in \{8,12\}$ leads to an equal number of iterations. It is also worth noting that larger n such values of the degree k have little effect on the overall computational complexity.

7. Concluding Remarks

We analyzed the non-overlapping DD method in the general framework of problems which involve coupling at interfaces that are manifolds of lower dimensions. This naturally gives rise to interface conditions formulated in fractional Sobolev spaces. Thus, the construction of the method is based on the spectral equivalence of the energy norm associated with the related Steklov-Poincaré operator and the Sobolev norm of index 1/2. The matrix formulation of the preconditioner uses the spectral equivalence of the Schur complement and the discrete fractional Laplacian over the interface $\gamma$, denoted by ${\mathbb{L}}_{\gamma }$. For a long time, the direct application of this approach was avoided due to the lack of efficient numerical methods for solving fractional diffusion problems in domains (interfaces) of general geometry. The situation has changed significantly in recent years [11].

The proposed new non-overlapping DD preconditioner is based on BURA approximation of ${\mathbb{L}}_{\gamma }^{1/2}$. A complete theoretical analysis of the method is provided, including, in particular, a characterization of the introduced scaling parameters ${\sigma }_1$ and ${\sigma }_2$. Finally, we have proved the overall optimality of the computational complexity. The presented numerical results illustrate in detail the theoretical estimates, ending with some practical tips for choosing the degree of rational approximation k and the parameter ${\sigma }_1{\sigma }_2$.

Although presented in 2D formulation, the theory can easily be extended to the 3D case. Note that the estimates of the condition number of $C_{0.5,k}^{BURA}\left({\mathbb{L}}_{\gamma }\right)$ are robust, with respect to the spatial dimension and the geometry of $\gamma$. This argues that we can get optimal preconditioning properties of the BURA based non-overlapping DD methods for a very wide class of problems. This statement is supported by the numerical results published in [1], where exact solution of the systems with fractional Laplacian over the interface $\gamma$ is used. Although this approach is computationally too expensive, the reported results are indicative for what we can expect in terms of BURA-based DD behavior in cases of general geometry of the interface and heterogeneous diffusion coefficients in 2D and 3D. We also note that a coarse grid solver may be needed to stabilize the DD iterations for a larger number of subdomains.

Possibilities for generalizing the presented results include: theoretical and experimental analysis of the new method for three-dimensional problems; development of robust domain decomposition preconditioners for the case of strongly heterogeneous and anisotropic media; computational complexity analysis for the case of fractal geometry interfaces. The application of the BURA-based non-overlapping domain decomposition approach to couple non-linear problems is a separate challenging direction for obtaining future results with great potential and practical significance.

One of the remarkable recent examples of effective application of domain decomposition approaches to elliptic PDEs is that of the Poisson-Boltzmann equations (see [16,17]) widely arising in the colloid science, protein electrostatics, implicit solvent models, and plasma physics. Development of BURA based non-overlapping domain decomposition preconditioners for the Poisson-Boltzmann elliptic boundary value problems gives another good perspective for state-of-the-art real-world applications.