How to Best Choose the Outer Coarse Mesh in the Domain Decomposition Method of Bank and Jimack

In Ciaramella et al. (2020) we defined a new partition of unity for the Bank–Jimack domain decomposition method in 1D and proved that with the new partition of unity, the Bank–Jimack method is an optimal Schwarz method in 1D and thus converges in two iterations for two subdomains: it becomes a direct solver, and this independently of the outer coarse mesh one uses! In this paper, we show that the Bank–Jimack method in 2D is an optimized Schwarz method and its convergence behavior depends on the structure of the outer coarse mesh each subdomain is using. For an equally spaced coarse mesh its convergence behavior is not as good as the convergence behavior of optimized Schwarz. However, if a stretched coarse mesh is used, then the Bank–Jimack method becomes faster then optimized Schwarz with Robin or Ventcell transmission conditions. Our analysis leads to a conjecture stating that the convergence factor of the Bank–Jimack method with overlap L and m geometrically stretched outer coarse mesh cells is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$1-O(L^{\frac {1}{2 m}})$\end{document}1−O(L12m).


Introduction
In 2001, Randolph E. Bank and Peter K. Jimack [3] introduced a new domain decomposition solver for the Bank-Holst adaptive meshing paradigm [2] for the adaptive solution of elliptic partial differential equations. The novel feature of the Bank-Jimack method (BJM) is that each of the subproblems is defined over the entire domain, but outside of the subdomain, a coarse mesh is used. A variant of this new domain decomposition solver using an augmented Lagrange multiplier technique is analyzed in [5] in the context of the abstract Schwarz framework. The BJM in [3] is formulated as a residual correction method, and it is not easy to interpret how and what information is transmitted between subdomains through the outer coarse mesh each subdomain has. A similar difficulty of interpretation existed as well for Additive Schwarz and Restricted Additive Schwarz [11,15,25]. This is very different compared to classical domain decomposition methods where this is well understood: classical Schwarz methods [24] exchange information through Dirichlet transmission conditions and use overlap, FETI [12,13] and Neumann-Neumann methods [6,22,23] use Dirichlet and Neumann conditions without overlap, and optimized Schwarz methods (OSMs), which go back to Lions, [20] use Robin or higher order transmission conditions and work with or without overlap, see [14,15] for an introduction and historic perspective of OSMs. In [10], we showed for a one-dimensional Poisson problem and two subdomains that if one introduces a more general partition of unity, then the BJM becomes an optimal Schwarz method, i.e. a direct solver for the problem converging in two iterations, and this independently of how coarse the outer mesh is. The BJM thus faithfully constructs a Robin type transmission condition involving the Dirichlet to Neumann map in 1D. We analyze here the BJM for the Poisson equation in 2 dimension and two subdomains, and show that with the modified partition of unity, the method can be interpreted as an OSM. Its convergence now depends on the structure of the outer coarse mesh each subdomain uses. In case of equally spaced coarse meshes, we prove that the asymptotic convergence factor is not as good as for an OSM. If one uses however a stretched coarse mesh, i.e. a mesh which becomes gradually more and more coarse in a specific way as one gets further away from the subdomain boundary, the method converges faster than the classical zeroth and second-order OSMs. Based on extensive numerical and asymptotic studies of the analytical convergence factor and the position of coarse points, we conjecture an asymptotic formula for the contraction factor of the BJM. Our analysis also indicates a close relation of the BJM to the class of sweeping type preconditioners [19], since the outer coarse mesh can be interpreted as an implementation of a PML transmission condition, but the BJM is not restricted to sequential decompositions without cross points.
Our paper is organized as follows: in Section 2, the BJM is recalled for a general PDE problem and its generalization by a partition of unity function is introduced (for the influence of partitions of unity on overlapping domain decomposition methods, see [16]). Moreover, for the Laplace problem in two dimensions the BJM is described in detail. The convergence analysis of the BJM is carried out in Section 3, where it is proved to be equivalent to an OSM. This important relation allows us to obtain sharp convergence results. Section 4 is devoted to extensive numerical experiments leading to our conjecture. Finally, our conclusions are presented in Section 5.

The Bank-Jimack Domain Decomposition Method
In this section, we give a precise description of the BJM, and introduce our model problem and the Fourier techniques that we will use.

Description of the Method
Let us consider a general self-adjoint 1 linear elliptic PDE Lu = f in a domain Ω with homogeneous Dirichlet boundary conditions on ∂Ω. Discretizing the problem on a global fine mesh leads to a linear system Au = f , where the matrix A is the discrete counterpart of L, u is the vector of unknown nodal values on the global fine mesh, and f is the load vector.
To describe the BJM, we decompose Ω into two overlapping subdomains, Ω = Ω 1 ∪Ω 2 . The unknown vector u is partitioned accordingly as u = u 1 , u s , u 2 , where u 1 is the vector of unknowns on the nodes in Ω 1 \ Ω 2 , u s is the vector of unknowns on the nodes in the overlap Ω 1 ∩ Ω 2 , and u 2 is the vector of unknowns on the nodes in Ω 2 \ Ω 1 . We can then write the linear system Au = f in block-matrix form, The idea of the BJM is to consider two further meshes on Ω, one identical to the original fine mesh in Ω 1 , but coarse on Ω \ Ω 1 , and one identical to the original fine mesh in Ω 2 , but coarse on Ω \ Ω 2 . This leads to the two further linear systems where we introduced the restriction matrices M j , j = 1, 2, to map the fine-mesh vectors f j to the corresponding coarse meshes, and I 1 and I 2 are identities of appropriate sizes. Notice that, depending on the chosen discretization scheme, one could get C 1 = B 1 and C 2 = B 2 , which leads to symmetric matrices A Ω 1 and A Ω 2 . However, this symmetry is not generally guaranteed, as we are going to see in the next sections. The BJM as a stationary iteration is then described by Algorithm 1.
In [10] we studied the BJM for a one-dimensional problem and showed that, in general, it does not lead to a convergent stationary iteration. 2 To correct this behavior we introduced a discrete partition of unity D 1 + D 2 = I , where I is the identity matrix and D 1 and D 2 are two matrices that for a one-dimensional problem must have the form (× denote arbitrary entries satisfying the sum condition) Using these matrices, we modified the BJM by replacing Step 2.3 in Algorithm 1 with ⎡ ⎣ u k+1 This leads to an iterative method that we proved to be convergent and equivalent to an optimal Schwarz method [18] for the one-dimensional Poisson problem [10]. In [10] we also showed, by direct numerical experiments, that this equivalence does not hold for the two-dimensional Poisson problem. Our goal here is to analyze the convergence of the BJM for two-dimensional problems. Notice that, in what follows, we always refer to BJM as the method obtained by using (2.4) in Algorithm 1.

The BJM for the Poisson Equation in 2D
Let us consider the problem where is the Laplace operator and f is a sufficiently regular right-hand side function. We consider a uniform grid in Ω = (0, 1) × (0, 1) of N interior points in each direction and , and two partially coarse meshes corresponding to the left subdomain Ω 1 (middle row) and to the right subdomain Ω 2 (bottom row). The black dots represent the x-coordinates of the interfaces, namely mh and h mesh size h := 1 N+1 ; see, e.g., Fig. 1a. We then discretize (2.5) by a second-order finitedifference scheme, which leads to a linear system Au = f , where A ∈ R N 2 ×N 2 is the classical (pentadiagonal) discrete Laplace operator. This system can be easily partitioned as in (2.1). We assume that the vector of unknowns u is obtained as u = [u 1 , . . . , u N ] , where u j ∈ R N contains the unknown values on the j th column of the grid. In this case, the matrix A can be expressed in the Kronecker format A = I y ⊗ A x + A y ⊗ I x , where A x and A y are N × N one-dimensional discrete Laplace matrices in directions x and y, and I x and I y are N × N identity matrices.
The BJM requires two partially-coarse grids. We assume that, in direction x our decomposition Ω = Ω 1 ∪ Ω 2 has n 1 interior points Ω 1 \ Ω 2 , n 2 interior points in Ω 2 \ Ω 1 , and n s points in Ω 1 ∩ Ω 2 ; see Fig. 1a. For our analysis, the coarsening is performed only in x-direction, 3 as shown in Fig. 1b, while the grid in direction y is maintained fine. The partially coarse-grids have m 2 coarse points in Ω 2 \ Ω 1 (Fig. 1b, middle row) and m 1 coarse points in Ω 1 \ Ω 2 (Fig. 1b, bottom row) and the corresponding mesh sizes are h 1 and h 2 . 4 If we denote by A x,1 and A x,2 the one-dimensional finite-difference Laplace matrices in x-direction, then the partially-coarse matrices A Ω 1 and A Ω 2 are where I x,1 and I x,2 are identities of sizes n 1 + n s + m 2 and m 1 + n s + n 2 , respectively. Notice that, the matrices A x,1 and A x,2 are classical second-order finite difference matrices in 1D, defined on the union of two uniform grids. Therefore, the only entries that differ from a standard finite-difference formula are the ones corresponding to the stencil across the mesh changes. For example, the five-point formulas for A x,1 on fine and coarse meshes , for j ≥ n 1 + n 2 + 2, while at the point across the mesh change (see also Fig. 1), we have for j = n 1 + n s + 1.
The matrices A Ω 1 and A Ω 2 can be partitioned exactly as in (2.2), and the restriction matrices T 1 and T 2 have now the forms where M 1 ∈ R m 1 ×n 1 and M 2 ∈ R m 2 ×n 2 are one-dimensional restriction matrices and I n s +n 2 and I n 1 +n s are identity matrices of sizes n s + n 2 and n 1 + n s . It remains to describe the matrices D 1 ∈ R Nn s ×Nn s and D 2 ∈ R Nn s ×Nn s used in (2.4). These form a partition of unity, that is D 1 + D 2 = I Nn s , where I Nn s is an identity of size Nn s , and have the forms We have then described all the components that allow us to use the BJM (namely Algorithm 1) for the two-dimensional Poisson problem (2.5). The choice of discretization by the finite-difference method allows us to perform a detailed convergence analysis based on the diagonalization obtained in Section 2.3.

A Discrete Fourier Expansion and the η − Equation in 1D
The finite-difference matrices A, A Ω 1 and A Ω 2 have similar structures based on Kroneckerproduct expansions: the matrix components in direction y are the same and are not coarsened. Hence, the one-dimensional discrete Laplace matrix A y appears unchanged in A, A Ω 1 and A Ω 2 , while the matrix A x appearing in A is replaced in A Ω 1 by A x,1 and in It is important to notice that A y is a tridiagonal Toeplitz matrix having values 2/h 2 on the main diagonal and values −1/h 2 on the first upper and lower diagonals. It is well-known that A y can be diagonalized as U A y U = , where = diag(λ 1 , . . . , λ N ) with λ j > 0, and the columns of the orthogonal matrix U ∈ R N×N are normalized discrete Fourier sine modes. If one now defines U := U ⊗ I x , then it is possible to block-diagonalize A, where we used the property ( for any matrices C 1 , C 2 , C 3 , and C 4 such that the matrix products C 1 C 3 and C 2 C 4 can be formed. This is the discrete version of a Fourier sine diagonalization of the continuous problem (2.5); see, e.g, [8]. Notice that each component u j ∈ R N still represents a vector of nodal values on the j th row of the discretization grid. Hence, we can decompose it as u j = [ u j,1 , u j,s , u j,2 ] , where u j,1 ∈ R n 1 has values on the nodes in Ω 1 \ Ω 2 , u j,2 ∈ R n 2 has values on the nodes in Ω 2 \ Ω 1 , and u j,s ∈ R n s has values on the nodes in Ω 1 ∩ Ω 2 . Now, using the block-diagonalized form (2.9)-(2.10), we will rewrite the BJM algorithm for each component j of u. Given an approximation u k obtained at the kth iteration of Algorithm 1, one can compute u k = U u k and r k = U r k and rewrite Step 2.1 as Similarly as for the system Au = f , we can transform the residual subsystems of Step 2.2.
To do so, we define U i := U ⊗ I x,i , for i = 1, 2, such that v k = U 1 v k and w k = U 2 w k , and write the subsystems A Ω 1 v k = T 2 r k and A Ω 2 w k = T 1 r k as which allows us to obtain Now, using the structures of A Ω i given in (2.6), we obtain for i = 1, 2, and recalling the matrices T i , defined in (2.7), we get and Replacing (2.13), (2.14) and (2.15) into (2.12), we rewrite the residual systems in Step 2.2 as where D e 1 ∈ R N×(n 1 +n s +m 2 ) and D e 2 ∈ R N×(m 1 +n s +n 2 ) are given by Now, using (2.17) we get where u j is the j th coefficient of the Fourier sine expansion of u, f j is the j th Fourier coefficient of f , and η j = (πj ) 2 . Hence, if we would know a continuous representation of the BJM for the solution to (2.19), then we could perform a Fourier convergence analysis similarly as it is often done at the continuous level for other one-level domain decomposition methods; see, e.g., [7,8,14]. This is exactly the focus of Section 3, where we will show that the BJM for the one-dimensional η − boundary value problem is an OSM. This equivalence will allow us to perform a detailed Fourier convergence analysis of the BJM.

Convergence Analysis of the BJM
Motivated by the results in Section 2.3, we study now the BJM for the solution of a onedimensional discrete η − problem and prove that this is equivalent to a discrete OSM, see for example [26]. Our analysis will reveal that the BJM produces implicitly some particular Robin parameters, dependent on η, in the equivalent OSM. Since the chosen discretization for the OSM is consistent and convergent, one can pass to the limit from the discrete to the continuous level. Therefore, we will obtain that the continuous limit of the BJM is an OSM, where the Robin parameters are the continuous limits of the discrete Robin parameters of the BJM. Once this equivalence interpretation is established, we will study the dependence of the continuous convergence factor of the BJM with respect to η (hence the Fourier frequency), to the size of the overlap, to the number of coarse points and their location. The main steps of the described analysis are organized in four subsections. In Section 3.1 we recall the OSM, derive its convergence factor at the continuous level and then obtain a discretization based on the finite-difference method for non-uniform grids; see, e.g., [27]. In Section 3.2, we show the equivalence between the BJM and the discrete OSM and discuss the BJM convergence factor in the continuous limit. Sections 3.3 and 3.4 focus on the analysis of the BJM convergence factor for uniform and non-uniform coarse grids.

The OSM for the One-Dimensional η − Equation
To recall the OSM for we consider an overlapping domain decomposition (0, 1) = (0, β := x ) ∪ (α := x m , 1); see Fig. 1b, top row. Given an appropriate initialization pair (u 0 1 , u 0 2 ), the OSM for (3.1) is 2) for k = 1, 2, . . . , where p 12 and p 21 are two positive parameters that can be optimized to improve the convergence of the iteration; see, e.g., [14]. This optimization process gives the name Optimized Schwarz Method to the scheme (3.2). In fact, the convergence factor of the method depends heavily on p 12 and p 21 . To compute this convergence factor, we can assume that f = 0 (working by linearity on the error equations). The general solution of the first subproblem in (3.2) , and since u 2 (1) = 0, we find that ). Using the Robin transmission condition at x = α in the second subproblem of (3.2), we find which leads to Similarly, using the Robin condition at the point x = β in the first subproblem of (3.2) we find shows that the convergence factor over a double iteration of the OSM is Notice that the convergence factor ρ depends on η, the two Robin parameters p 12 and p 21 , and on the positions of the interfaces α and β (hence the length of the overlap L := β − α).
To obtain a discrete formulation of the OSM, we consider two uniform grids of size h in the subdomains (0, β) and (α, 1) as the ones shown in Fig. 1b, top row. Using the finitedifference method applied to these grids, we discretize the two subproblems in (3.2) and obtain the linear systems where A OSM,j ∈ R (n j +n s )×(n j +n s ) and f j ∈ R n j +n s , j = 1, 2, are and the matrices F 1 ∈ R n 1 +n s ×n 2 +n s and F 2 ∈ R n 2 +n s ×n 1 +n s are such that for any g ∈ R n 2 +n s and h ∈ R n 1 +n s . Notice that, since η, p 12 , p 21 > 0 for any h > 0 the matrices A OSM,1 and A OSM,2 are strictly diagonally dominant, hence invertible. Therefore, the OSM (3.6) is a stationary method whose standard form (see, e.g., [9]) is . This is sometimes also called an optimized block Jacobi algorithm; see, e.g., [26]. If convergent, this iterative procedure generates a sequence that converges to the solution of the augmented problem In our analysis, another equivalent form of the the discrete OSM (3.7) will play a crucial role. This is the so-called optimized restricted additive Schwarz (ORAS) method, which is defined as where r k = f − A u k , R 1 ∈ R (n 1 +n s )×N , and R 2 ∈ R (n 2 +n s )×N are restriction matrices of the form while R 1 ∈ R (n 1 +n s )×N and R 2 ∈ R (n 2 +n s )×N are similar restriction matrices, but corresponding to a non-overlapping decomposition satisfying R 1 R 1 + R 2 R 2 = I N ; see [26] for more details. It is proved in [26] that (3.8) and (3.7) are equivalent for any R 1 and R 2 , as the ones considered in this section, that induce a consistent matrix splitting.

The BJM as an OSM for the One-Dimensional η − Equation
Let us first recall the BJM for the one-dimensional problem (3.1) and state explicitly all the matrices that we need for our analysis. We consider the grids shown in Fig. 1b and the finitedifference method for non-uniform grids; see, e.g., [27]. The full problem on the global fine mesh (Fig. 1b, top row) is where A ∈ R N×N is a tridiagonal symmetric matrix that we decompose as The matrices A 1 ∈ R n 1 ×n 1 , A s ∈ R n s ×n s , and A 2 ∈ R n 2 ×n 2 , are tridiagonal and have the form while B 1 ∈ R n 1 ×n s and B 2 ∈ R n 2 ×n s are zero except for one corner entry: Hence for a given approximation u k = [(u k 1 ) , (u k s ) , (u k s ) ] , the residual r k is The correction problems on the two partially coarse grids (Fig. 1b, middle and bottom rows), are where A Ω 1 ∈ R (n 1 +n s +m 2 )×(n 1 +n s +m 2 ) , A Ω 2 ∈ R (n 2 +n s +m 1 )×(n 2 +n s +m 1 ) , T 1 ∈ R (n 2 +n s +m 1 )×N , and T 2 ∈ R (n 1 +n s +m 2 )×N have the forms given in (2.2), with A 1 , A s , A 2 , B 1 , and B 2 as above. The matrices C 1 ∈ R n s ×m 1 and C 2 ∈ R n s ×m 2 are The matrices A 1 ∈ R m 1 ×m 1 and B 1 ∈ R m 1 ×n s in the BJM method in Algorithm 1 are We do not need to specify the restriction matrices M 1 and M 2 , because they multiply the residual components r 1 and r 2 , which are zero as shown in the upcoming Lemma 3.1. The matrices M j do not play any role in the convergence of the method if our new partition of unity is used. However, if the original partition of unity proposed in [3] is considered, then they contributes to the convergence behavior. Finally, the partition of unity diagonal matrices D 1 ∈ R n s ×n s and D 2 ∈ R n s ×n s have the structures given in (2.3). Notice that, since η > 0, the tridiagonal matrices A Ω 1 and A Ω 1 are strictly diagonally dominant for any The BJM in Algorithm 1 consists of iteratively computing the residual (3.11), solving the two correction problems (3.12) and then computing the new approximation using (2.4).
We are now ready to prove the equivalence between the BJM and the discrete OSM. To do so, we need an important property of the BJM proved in the next lemma.

Proof
We only sketch the proof here, since the result is proved in detail in [10]. Moreover, we consider only r k 1 , because the proof for r k 2 is similar. Using equations (3.11) and (2.4), we compute because of (3.12) at k − 1. Now using the structures of B 1 , D 1 and D 2 we get independently of the middle elements of D 1 , 5 and thus By a similarly calculation, one can show that B 1 D 2 w k−1 s = 0, also independently of the middle elements of D 2 , which proves that r k 1 = 0 for k = 1, 2, . . .
are well-defined and we can compute the entries we need for our analysis using the following lemma. 6 Lemma 3.2 The first element of the inverse of the n × n tridiagonal matrix Proof The first element of the inverse of T is the first component u 1 of the solution of the linear system The solution satisfies the recurrence relation −u j +1 + au j − u j −1 = 0, j = 2, 3, . . . , n− 1, whose general solution is u j = C 1 λ j 1 + C 2 λ j 2 with λ 1,2 the characteristic roots of λ 2 − aλ + 1 = 0 given in the statement of the lemma. The two boundary conditions to determine the constants C 1,2 are 2 ) + a(C 1 λ n 1 + C 2 λ n 2 ) = 0. Solving this linear system for C 1,2 gives (using that 3−i = 2 if i = 1 and 3−i = 1 if i = 2) Inserting these constants into u j and evaluating at j = 1 gives , which upon simplification, using the Vieta relations satisfied by the roots, i.e. λ 1 λ 2 = 1 and λ 1 + λ 2 = a, leads to the result.

Lemma 3.3
The matrices C 2 A −1 2 B 2 and C 1 A −1 1 B 1 are given by with the function μ(n) from Lemma 3.2.
Proof For the first result, using the sparsity patterns of C 2 and B 2 , we obtain and we thus need to find the first entry of A −1 2 . Defining a 1 := 2h 1 h + ηh 2 1 , b 1 := −2h 1 h+h 1 , and a := 2 + ηh 2 1 , and multiplying by h 2 1 , we obtain precisely a matrix like in Lemma 3.2, and therefore ((h 2 1 A 2 ) −1 ) 11 = μ(m 2 ), which shows the first claim. For the second one, it suffices to notice that Lemma 3.2 also holds if the matrix is reordered from top left to bottom right, and can thus be used again.
Notice that the matrices C 2 A −1 2 B 2 and C 1 A −1 1 B 1 are operators which are non-zero only on the last column of the overlap. Therefore, as we are going to see in the rest of this section, they can be interpreted as generalized Robin boundary conditions. Now, using the Schur- , we can introduce the matrices A 1 and A 2 : which allow us to prove the following result.

Lemma 3.4 The matrices A 1 and A 2 are invertible and the inverses of A Ω 1 and A Ω 2 have the forms
and

14)
where Proof We prove the result for A 1 . The proof for A 2 can be obtained exactly by the same arguments. Recalling that η > 0, a direct inspection of the matrix A Ω 1 reveals that it is strictly diagonally dominant. Hence, det(A Ω 1 ) = 0. Now, consider the block structure of A Ω 1 given in (2.2). Since A 2 is invertible, we factorize A Ω 1 as ⎡ where I n 1 , I n s , and I m 2 are identity matrices of sizes n 1 , n s , and m 2 . This factorization allows us to write 0 = det(A Ω 1 ) = det( A 2 )det( A 1 ), which implies that det( A 1 ) = 0. Now, a straightforward calculation using the previous factorization allows us to get (3.13).
Now, we notice that the BJM can be written (using (3.12) and (2.4)) in the compact form where the block-diagonal matrices T 1 ∈ R (n 1 +n s +m 2 )×N and T 2 ∈ R (m 1 +n s +n 2 )×N are A direct calculation using Lemma 3.1 (hence that r k 1 = 0 and r k 2 = 0) and Lemma 3.4 (hence the formulas (3.13) and (3.14)) allows us to obtain where the matrices R 1 and R 2 are the ones given in (3.9). Since the results proved in Lemma 3.1 are independent of the middle diagonal entries of D 1 and D 2 , we can choose them such that the equalities are fulfilled. Therefore, the BJM (3.15) becomes Similarly, A 2 and A OSM,2 are equal except for the top-left corner elements, which are Therefore, if one chooses (3.18) then A j = A OSM,j for j = 1, 2. Replacing this equality into (3.17), we obtain that the BJM is equivalent to the ORAS method (3.8), and hence to the discrete OSM (3.6). We summarize our findings in the following theorem. Notice that Theorem 3.5 has the following important consequence. Since the discrete OSM (3.6) is obtained by a consistent and convergent discretization of the continuous OSM (3.2), we find that, in the limit for h → 0, the continuous counterpart of the BJM is the OSM (3.2). This will allow us to study in Sections 3.3 and 3.4 the convergence factor of the BJM at the continuous level. For this purpose, from now on, we denote by p 12 (h, η, h 1 ) and p 21 (h, η, h 2 ) the two Robin parameters of (3.18) to stress their dependence on the discretization size h, the (Fourier) parameter η and the coarse mesh sizes h 1 and h 2 . Notice that μ(m 2 ) and μ(m 1 ) in (3.18) depend on h, h 1 , h 2 and η (see Lemma 3.3). Recalling the results obtained in Section 3.1, the continuous BJM convergence factor is given by (3.5), where p 12 and p 21 are the limits for h → 0 (with m 1 and m 2 fixed) of the parameters chosen in Theorem 3.5.
It is important to remark at this point that the first coarse points, namely the point n 1 + n s +1 for the first mesh and the point m 1 for the second mesh, are located at distance h from the interfaces. With this choice we were able to define discrete finite-difference derivatives across these points and in Sections 3.3 and 3.4 we will take limits for h → 0, while keeping the numbers m 1 and m 2 of the coarse points fixed.
Finally, we wish to remark that all the calculations performed in this section, except for the precise formulas for μ(m 2 ) and μ(m 1 ) in Lemma 3.3, remain valid if, instead of uniform coarse grids, one considers two coarse grids which are non-uniform, in the sense that the m 1 points in Ω 1 \ Ω 2 and the m 2 points in Ω 2 \ Ω 1 are not uniformly distributed, leading to invertible matrices A 1 and A 2 . Therefore, the equivalence between BJM and OSM remains valid also in the case of non-uniform coarse grids.

Uniform Coarse Grid
The goal of this section is to study the contribution of uniform coarse grids to the convergence of the BJM for the solution to (2.5). For simplicity, we assume that the two partially coarse grids have the same number of coarse points m := m 1 = m 2 . To satisfy this condition, we fix the size of the overlap L and choose α = 1−L 2 and β = 1+L 2 . In this case, we also have that h 1 = h 2 = 1−β−h m . We consider the cases of m = 2, m = 3, and m = 4 coarse points.
For the sake of clarity, we first summarize the structure of our analysis. For each given m ∈ {2, 3, 4}, we first consider the corresponding BJM Robin parameters, whose explicit formulas can be obtained as in Lemma 3.3, and then pass to the limit for h → 0 to get their continuous counterparts. These continuous parameters will be replaced into the formula (3.5), which will give us the continuous convergence factor of the BJM corresponding to the given m, to a fixed (Fourier) parameter η, and to the size of the overlap L. For fixed m and given values of L we will numerically compute the maximum of the convergence factor with respect to the (Fourier) parameter η. This will allow us to study the deterioration of the contraction factor for decreasing size L of the overlap. While performing this analysis, we compare the convergence of the BJM to the one of the OSM with optimized parameter.
From the convergence factor ρ of the OSM in (3.5), we see that choosing gives ρ = 0 for the frequency η. These are thus the optimal parameters for this frequency, and make the OSM a direct solver for the corresponding error component. For m = 2 coarse points, proceeding as in the proof of Lemma 3.2 to compute the corresponding μ(m 2 ) = μ(m 1 ) and using (3.18), we get the (discrete) BJM Robin parameters  Recalling that h 1 = h 2 = 1−β−h 2 and taking the limit for h → 0, we obtain (3.20) We see that the Robin parametersp 12 andp 21 are rational functions of the Fourier parameter η with coefficients depending on the outer subdomain sizes 1 − β and α. In Fig. 2, we compare the Robin parameterp 12 of the BJM for m = 2 (blue line) with the optimal Robin parameter p * 12 of the OSM (black dashed line) for three different values of the overlap L. We observe that for small η the Robin parameters of both methods are quite close, which indicates that the BJM method performs well for low-frequency error components. This is clearly visible in Fig. 3, where we plot the corresponding convergence factors (as functions of η) insertingp 12 andp 12 into (3.5) 7 for two different overlaps L, using α = 1−L 2 and β = 1+L 2 . We also see that the convergence factor clearly has a maximum at some η 2 (L), whose corresponding error mode converges most slowly, and convergence deteriorates when L becomes small. In Fig. 4 (left), we present the value η 2 (L) as functions of L and observe that it grows like O(L −1 ). The corresponding contraction factor, namelyρ 2 (L) := max η ρ 2 (η, L) := max η ρ η,p 12 (η, L),p 21 (η, L), α = 1−L 2 , β = 1+L 2 , is shown as function of L in Fig. 4 (right-dashed blue line, represented as 1 − ρ 2 (L)). Here, one can observe clearly that as L gets smaller the convergence deteriorates with an order O(L 1/2 ).
Let us now discuss the behaviorρ 2 (L) = 1 − O(L 1 2 ) shown in Fig. 4 (right): it was proved in [14] that the convergence factor of the OSM with overlap L behaves like ρ OSM = ) with second-order (Ventcell) transmission conditions. Hence, the OSM performs better than the BJM with a 7 To evaluate the convergence factor numerically one needs to factor out an exponential in the hyperbolic trigonometric functions to avoid overflow. uniform coarse grid with m = 2 uniformly distributed coarse points, 8 since convergence deteriorates more slowly when the overlap L goes to zero.
We have seen that, for only two points the BJM is already a good method for low frequencies, since the parametersp 12 andp 21 are very close to the optimal ones p * 12 and p * 21 for relatively small η. However, the convergence factor deteriorates with L faster than for the OSM. It is natural to ask: does the behavior of the BJM improve if more coarse points are used? The answer is surprisingly negative! In fact, the convergence factor remains of . To see this, we now repeat the analysis for uniform coarse grids with m = 3 and m = 4 points. For m = 3, we find the analog of (3.19) with E 2 (h) replaced by and for the corresponding optimized parameters when h goes to zero the analog of (3.20) with the rational function R 2 (L) replaced by In Fig. 2 (red lines) we show the Robin parameters of the BJM with m = 3 coarse points as a function of η and we compare it to the optimal Robin parameters of the OSM. We observe that they are closer compared to the m = 2 point case. This seems to suggest an improvement of the convergence factor, but the plots of the convergence factor in Fig. 3 show that this improvement is only minor compared to the case of m = 2 coarse mesh points. This is also confirmed by the results in Fig. 4  We thus conclude that the convergence factor of the BJM with a uniform coarse grid always behaves as 1 − O(L 1 2 ) independently of the number of coarse points of the grids. This shows that the OSM has a better convergence factor compared to the BJM with uniform coarse grids since its convergence factor behaves as 1 − O(L 1 3 ), but BJM with uniform coarse grids converges better than classical Schwarz, which has a convergence factor 1 − O(L), see [14]. Is the uniformity of the coarse grids the limiting factor for BJM? We address this in the next section.

Stretched Coarse Grid
We now consider stretched coarse grids, and start with m = 2 non-uniformly distributed coarse points with grid sizes h 1 1 , h 2 1 , h 1 2 , and h 2 2 , see Fig. 5 (second and third rows). Using the finite-difference method, we discretize our problem and obtain the two linear systems  1 and A Ω 2 have the block-structures given in (2.2) with the blocks corresponding to the coarse parts of the grids that are Proceeding as in Section 3.3 we find after some calculations discrete BJM parameters of the form (3.19), but with E 2 (h) replaced bỹ We now use the relations h 2 which shows that the coefficients in the rational function in η can now be controlled by the mesh parameterh 1 ! To understand the impact of this new degree of freedom from the coarse mesh, we assume for simplicity that α = 1−L 2 and β = 1+L 2 , and h 1 1 = h 1 2 and h 2 1 = h 2 2 . Insertingp 12 andp 21 into (3.5) and minimizing the maximum of the resulting convergence factor (3.5) over all frequencies η (using the MATLAB function fminunc), we find the best choice for the mesh stretching h 1 1 (L) that makes the convergence factor as small as possible. We show in Fig. 6 the behavior of the Robin parameterp 12 (η) (blue lines) compared to the OSM parameter p 12 (η) (black dashed lines) for different overlaps L. Clearly, the curves are very different from the ones corresponding to the uniform mesh (Fig. 2) which are very stable with respect to the overlap L. In the stretched case, the coarse mesh is strongly influenced by the overlap: the smaller the overlap, the more work needs/can be done in the optimization of the coarse points. The corresponding convergence factors are shown in Fig. 7 (blue lines), where one can now observe how they have two maxima. Hence, the optimization of the coarse points is solved when an equioscillation of the two maxima is obtained. If one compares these plots to the ones presented in Fig. 3, the enormous improvement obtained by optimizing the position of the m = 2 coarse points is clearly visible. This behavior is even more evident if one compares the deterioration ofρ 2 of Fig. 4 (right) with the corresponding one of Fig. 8   Finally, in Fig. 9 (left) we show the dependence of the frequencies η 1 and η 2 (the maximum points) on L and we observe that We prove these numerical observations in the next theorem. To solve this non-linear system asymptotically, we insert the ansatz h 1 1 = h 1 2 := C h 1 √ L and η 1 := C η 1 L − 1 2 and η 2 := C η 2 L − 3 2 into the system (3.22), expand for overlap L small and find the relations The solution is C η 1 = C η 2 = 8 and C h 1 = 1 2 , which leads when inserted with the ansatz into ρ 2 (η 1 , L) to (3.21) after a further expansion for L small.
We thus conclude that the convergence factor of the BJM with an optimized stretched coarse mesh with m = 2 points behaves better than the convergence factor of the OSM with Robin transmission conditions which is ρ OSM = 1 − O(L . Similar calculations as before (see also [21]) lead after expanding for h going to zero to the continuous Robin parameters of the BJM (3.20) with the rational function R 2 (L) replaced bỹ We thus have now two parameters from the stretched mesh from each side to optimize the convergence factor! We set again α = 1−L 2 and β = 1+L 2 , and h j 1 = h j 2 , j = 1, 2, 3, and insertingp 12 andp 21 into the convergence factor (3.5) and minimizing the maximum of the resulting convergence factor over all frequencies η, we find the best choice for the mesh stretching h 1 1 (L), h 2 1 (L) that makes the convergence factor as small as possible, shown in Fig. 7 for a typical example in red. We notice that now three local maxima are present and equioscillate. In Fig. 8 (left), we show how the optimized choice of the stretched mesh parameters h 1 1 (L), h 2 1 (L) decay when the overlap L becomes small, and observe that Similarly, in Fig. 9 (middle) we find for the maximum points η 1 , η 2 , and η 3 the asymptotic behavior  , and η 1 := C η 1 L − 1 3 , η 2 := C η 2 L − 3 3 , η 2 := C η 2 L − 5 3 into the system (3.24), we can solve the system asymptotically for the constants when the overlap L becomes small, which leads to (3.23).
The analysis for m = 4 stretched coarse points follows the same lines, and we find after a longer computation for the continuous Robin parameters of the BJM (3.20) with the rational function R 2 (L) replaced by (see also [21] for details) with the numerator and denominator given bỹ Proof We proceed as in the proof of Theorems 3.6 and 3.7.
These results for optimized stretched coarse grids with m = 2, m = 3, and m = 4 points lead us to formulate the following conjecture: This result shows that one should choose a geometric coarsening related to the overlap to form the outer coarse grid leading to the best performance for the Bank-Jimack domain decomposition algorithm. A practical approach is to just take a geometrically stretched grid with respect to the overlap size, and then to sum the step sizes h j and scale the result to the size of the outer remaining domain, sayL, to get the actual mesh sizesh j to use, This direct geometric stretching including the last grid cell is preasymptotically even a bit better, as one can see in Fig. 10.

Numerical Experiments
In this section, we present numerical experiments to illustrate our theoretical results. We start with experiments for equally spaced coarse meshes, and compare their performance with the optimized geometrically stretched ones. We consider both a case of constant overlap L and a case where the overlap is proportional to the mesh size. We then also explore numerically the influence of coarsening the meshes in the direction tangential to the interface. In all these cases, we study the performance of the BJM as a stationary method and as a preconditioner for GMRES. We discretize the Poisson equation (2.5) (defined on a unit square Ω = (0, 1) 2 ) using n 2 (interior) mesh points where n = 2 − 1, for = 5, 6, 7, is the number of interior points on the global fine mesh in each direction (Fig. 1). The results corresponding to a uniform coarsening in direction x are presented in Section 4.1. Section 4.2 focuses on optimized stretched coarsening in direction x. Finally, in Section 4.3 we study the effect of the coarsening in both directions x and y.

Uniform Coarsening in Direction x
We start with the equally spaced coarse mesh case, coarsened only along the x axis. At first, we consider the case with a constant overlap L = 1 16 , which corresponds to n s = 3, 5, 9 for = 5, 6, 7, respectively. Moreover, to test the methods in the cases studied by our theoretical analysis, we consider m = 2, 3, 4 coarse mesh points. The results of the The former shows the decay of the error corresponding to the BJM as a stationary iteration, while the latter presents the decay of the GMRES residuals along the iterations. All the plots show that the effect of the number of coarse points on the convergence is very mild. This corresponds to the results discussed in Section 3.3 and shown in Fig. 4 (right): if the overlap L is constant, the contraction factor does not improve significantly if more (uniformly distributed) coarse points are considered. The same effect can be observed in the GMRES convergence. Now, we wish to study the effect of the new partition of unity proposed in [10] and constructed using (2.3). This was used in all the experiments discussed above. If we use the original partition of unity, we already know from [10] that the BJM does not converge as a stationary method. Therefore, we use it only as a preconditioner for GMRES and obtain the results depicted in Fig. 13. By comparing the results of this figure with the ones of Fig. 12, we see that the effect of the new partition of unity is tremendous: GMRES converges much faster and is very robust against mesh refinements. For further information on the influence of the partition of unity on Schwarz methods, see [16]. Now, let us now consider an overlap proportional to the mesh size, namely L = 2h, and repeat the experiments already described. The corresponding results are shown in Figs. 14, 15 and 16. As before, we observe that the BJM method (as stationary iteration and as

Stretched Coarsening in Direction x
In this section, we repeat the experiments presented in Section 4.1, but we optimize the position of the coarse mesh points by minimizing numerically the contraction factor (as in Section 3.4). We begin with the case of constant overlap L = 1 16 . The corresponding numerical results are shown in Figs. 17 and 18. These results show that optimizing the coarse mesh leads to a faster method which is robust against the mesh refinement. However, due to the constant overlap, there is only little improvement with respect to the constant coarsening case. To better appreciate the effect of the mesh optimization, we consider the case with overlap L = 2h. The corresponding results are shown in Figs. 19 and 20. By comparing these results with the ones of Figs. 14 and 15, one can see clearly the improvement of the BJM convergence: the number of iterations (for both stationary and preconditioned GMRES methods) are essentially halved in the case of finer meshes.

Coarsening in Direction x and y
We conclude our numerical experiments by studying the effect of a (uniform) coarsening in both x and y directions. As before, we consider both cases L = 1 16 and L = 2h. The results shown in Figs. 21, 22, 23 and 24 indicate that a coarsening in direction y does not have a significant impact on the convergence of the BJM method.

Conclusions
We provided a detailed convergence analysis of the Bank-Jimack domain decomposition method for the Laplace problem and two subdomains. Our analysis reveals that one should coarsen the outer mesh each subdomain uses in a geometric progression related to the size of the overlap if one wants to get good convergence, and arbitrarily weak dependence on the overlap size is possible (see also [17] for a different technique reaching this). In order for these results to hold one has to use a slightly modified partition of unity in the Bank-Jimack algorithm, without which the convergence of the method is much worse. We obtained our results by an asymptotic process as the subdomain mesh size goes to zero, and thus the results hold at the continuous level.
A possibility for further optimization at the discrete level is the observation that the maxima in the optimized method, shown in Fig. 7, occur for very high values of η which represent a Fourier frequency, and thus may lie outside of the frequencies representable on the mesh used. This can be seen quantitatively for example from the stretched case for m = 4, where the largest η 4 = O(L − 7 4 ), and the highest Fourier frequency can be estimated as η = O(h −1 ), see [14]. Hence, if the overlap is of the order of the mesh size, L = h, η 4 would be already much larger than what the grid can represent, and we see in fact from (3.25) that only half the number of bumps would need to be taken in consideration for the optimization.
Funding Open access funding provided by University of Geneva.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.