remarkRemark hypothesisHypothesis claimClaim A robust two-level overlapping
preconditionerC. Ye, S. Fu, E. Chung, and J. Huang

# A robust two-level overlapping preconditioner for Darcy flow in high-
contrast media††thanks: Submitted to the editors DATE. Jizu Huang’s work is
supported by National Key Research and Development Project of China (No.
2023YFA1011705) and National Natural Science Foundation of China (No.
12131002). Eric Chung’s work is partially supported by the Hong Kong RGC
General Research Fund (Project numbers: 14305222 and 14304021). Shubin Fu’s
work is supported by start up funding of Eastern Institute for Advanced Study.

Changqing Ye Department of Mathematics, The Chinese University of Hong Kong,
Shatin, Hong Kong SAR<EMAIL_ADDRESS>Shubin Fu Eastern Institute
for Advanced Study, Ningbo, China<EMAIL_ADDRESS><EMAIL_ADDRESS>Eric T. Chung Department of Mathematics, The Chinese
University of Hong Kong, Shatin, Hong Kong SAR<EMAIL_ADDRESS>Jizu Huang LSEC, Academy of Mathematics and Systems Science, Chinese Academy
of Sciences, Beijing 100190, China, and School of Mathematical Sciences,
University of Chinese Academy of Sciences, Beijing 100049, China ().
<EMAIL_ADDRESS>
###### Abstract

In this article, a two-level overlapping domain decomposition preconditioner
is developed for solving linear algebraic systems obtained from simulating
Darcy flow in high-contrast media. Our preconditioner starts at a mixed finite
element method for discretizing the partial differential equation by Darcy’s
law with the no-flux boundary condition and is then followed by a velocity
elimination technique to yield a linear algebraic system with only unknowns of
pressure. Then, our main objective is to design a robust and efficient domain
decomposition preconditioner for this system, which is accomplished by
engineering a multiscale coarse space that is capable of characterizing high-
contrast features of the permeability field. A generalized eigenvalue problem
is solved in each non-overlapping coarse element in a communication-free
manner to form the global solver, which are accompanied by local solvers
originated from additive Schwarz methods but with a non-Galerkin
discretization to derive the two-level preconditioner. We provide a rigorous
analysis indicates that the condition number of the preconditioned system
could be bounded above with several assumptions. Extensive numerical
experiments with various types of three-dimensional high-contrast models are
exhibited. In particular, we study the robustness against the contrast of the
media as well as the influences of numbers of eigenfunctions, oversampling
sizes, and subdomain partitions on the efficiency of the proposed
preconditioner. Besides, strong and weak scalability performances are also
examined.

###### keywords:

preconditioner, domain decomposition method, heterogeneous Darcy flow,
multiscale space

65N55, 65F08, 65F10

##  1 Introduction

Simulating subsurface fluids through porous media is important for many
practical applications, such as reservoir simulation, nuclear water
management, and groundwater contamination predictions. In subsurface modeling,
properties of the medium such as permeability may vary by several orders of
magnitude spatially and have complicated structures. Moreover, high-resolution
models have become available due to advances in reservoir characterization
techniques. All these factors impose prohibitively high computational burdens
on existing simulators and have led to the developments of model reduction
techniques such as upscaling [50, 20, 6], in which people aim to homogenize
high resolution models to obtain reduced ones that can be solved with
permissible computational costs. Although upscaling has been successfully
applied for considerable problems, it usually fails to resolve all small-scale
information of the interested quantities. To overcome the bottlenecks of
upscaling, various types of multiscale methods including the multiscale finite
element method [33] with its extensions the mixed multiscale finite element
method [15, 1], Generalized Multiscale Finite Element Methods (GMsFEM) [21,
17], the multiscale finite volume method [39, 27] and the multiscale mortar
mixed finite element method [4] have been proposed. The core idea of these
methods is to solve target problems in coarse grids with carefully engineered
multiscale basis functions, which are usually solutions of some well-designed
sub-problems containing heterogeneity information of the original geological
models. Hence, those multiscale methods can usually capture small-scale
details with affordable computational resources. Nevertheless, the performance
of multiscale methods may deteriorate if the permeability fields exhibit
extremely high contrast and long channels [40], which necessitates developing
robust and efficient solvers for the original large-scale fine problems.

Direct solvers such as MUMPS [3, 2] apply a LU decomposition to the target
matrix and then solve the linear system efficiently for different right-hand
terms. However, in reservoir simulations, we need to update linear systems due
to the change of coefficient profiles, which implies the scenario of solving
multiple right-hand sides with one single matrix may not exist and the
advantage of direct solvers is greatly diminished. Besides, direct solvers
empirically require a huge amount of memory and perform poorly in terms of
scalability. Therefore, direct solvers are not good candidates for large-scale
flow simulations, which renders preconditioned iterative solvers a suitable
choice. Classical two-level Domain Decomposition (DD) preconditioners have
been implemented for various problems [47, 18]. In particular, the two-level
restricted additive Schwarz preconditioners combined with different solving
strategies can greatly improve and expand the scope of application of the
domain decomposition algorithm for large-scale flow problems [51, 53, 38, 52].
However, the iteration numbers of standard DD preconditioners theoretically
depend on the contrast ratio of the permeability field, as a result, the
performance for arbitrarily complicated media may not be guaranteed. One
remedy is to replace classic polynomial-based coarse spaces with non-standard
coarse spaces that embedded with heterogeneous models information [25, 46, 49,
42, 12, 34, 19, 35, 36, 41, 23, 24, 43, 44, 22, 37, 30, 31]. In [25, 49, 42],
multiscale basis functions [33] are utilized to construct two-level or two-
grid preconditioners for highly heterogeneous incompressible flow models, it
is also shown these preconditioners perform better than polynomial-based
classic preconditioners. However, numerical experience suggests that these
preconditioners are robust only if there are no long channels across coarse
grids. To deal with more challenging models, the idea of using spectral
functions was proposed and could be found in an extensive literature [12, 34,
19, 35, 36, 41, 23, 24, 43, 44, 22, 37, 30, 31]. It can now be confirmed from
these works that utilizing eigenfunctions of some carefully designed local
spectral problems is to some extent indispensable for constructing a robust
and efficient preconditioners. However, most of these works intend to solve
elliptic problems in second-order forms, few of them have been devoted to
studying the elliptic problems in first-order forms, via which we can generate
mass-conserving velocity fields that are valuable in long-time subsurface
simulations. Recently, a novel mixed GMsFEM [13, 28, 29] was proposed for
solving steady state incompressible flows in the first-order form. This method
contains several steps: (1) obtaining the diagonal velocity mass matrix using
the trapezoidal quadrature rule; (2) eliminating the velocity field to derive
a system with pressure as the only unknown; and (3) constructing multiscale
basis functions for pressure field as a model reduction. However, it is
observed in the numerical experiments, this method [13] requires using many
bases to achieve a high accuracy that is needed in solving the transport
equation for two-phase flow simulations.

An efficient and robust two-grid preconditioner within the framework of mixed
GMsFEM was proposed in [54] for solving single-phase and two-phase
incompressible heterogeneous flow problems. The robustness of this
preconditioner lies in the utilization of eigenfunctions of well-designed
spectral problems to form the coarse space. It is shown that the incorporation
of eigenfunctions in coarse spaces can tremendously improve performance over
standard two-grid preconditioners. However, the local spectral problems are
defined in strongly overlapped subdomains in [54], which greatly affects the
parallel efficiency.

In this paper, we develop a parallelized two-level overlapping preconditioner
that is based on GMsFEM for solving the fine-scale pressure system. In each
coarse element, a spectral problem is solved and eigenfunctions corresponding
to those small eigenvalues are selected to form a coarse space and then a two-
level preconditioner (the local solvers here are discretizations of zero-
Dirichlet problems defined on overlapped subdomains) is established for the
iterative solver. We emphasize that each local spectral problem is solved in a
coarse element that has no overlaps with other coarse elements, which saves
MPI inter-node communications and improves parallel efficiency compared with
previous efforts [12, 34, 19, 35, 36, 41, 23, 24, 43, 44, 22, 37, 30, 31, 54].
We here address preliminary results on robust DD preconditioners for the mixed
form of Darcy flow problems in [22], where the authors also demonstrate
extensions on Stokes’ and Brinkman’s equations. By introducing steam
functions, they obtain a positive definite bilinear form that is suitable for
constructions of eigenvalue problems. Our approach distinguishes from [22] in
that we sacrifice the pursuit of complete mixed form but focus on an
alternative positive semidefinite bilinear form w.r.t. pressure, which also
allows us to obtain a velocity field that satisfies mass conservation.
Additionally, our approach does not rely on partition functions, and the
eigenvalue problems are formulated on a non-overlapping decomposition,
resulting in a more straightforward implementation. Moreover, rich numerical
experiments for several different types of highly heterogeneous models with a
resolution up to $512^{3}$ are exhibited. In particular, we examine the
robustness of our preconditioner against contrast ratios of media. Strong and
weak scalability tests and performance comparisons with an out-of-the-box
preconditioner from PETSc [8] are provided. Besides, the influence of numbers
of eigenfunctions, oversampling sizes as well as subdomain partitions will be
also investigated. The proposed two-level preconditioner shows an encouraging
computational performance. Under several assumptions, we derive an a priori
bound of condition numbers of preconditioned linear systems. Our analysis is
featured for considering no-flux boundary conditions which result in singular
linear systems and novel designs of local solvers.

The rest of this paper is arranged as follows. Section 2 introduces the basic
model of subsurface flows, the definition of notations and fine-scale
discretization. Section 3 illustrates a systematic way to construct the coarse
space and the preconditioner. An analysis for the proposed preconditioner is
presented in Section 4. Numerical results are shown in Section 5. Finally, a
conclusion is given.

##  2 Preliminaries

###  2.1 The model problem and the mixed finite element discretization

We consider the following model problem for the unknown pressure field $p$ on
a bounded Lipschitz domain $\Omega$ in $\mathbb{R}^{d}$ where $d=2$ or $3$:

(1) $\left\\{\begin{aligned} -\operatorname{div}\left(\kappa\nabla
p\right)=&f\quad&&\text{in}\quad\Omega,\\\ \kappa\nabla
p\cdot\bm{n}=&0\quad&&\text{on}\quad\partial\Omega,\end{aligned}\right.$

where $\partial\Omega$ is a Lipschitz continuous boundary and $\bm{n}$ is the
unit outward normal vector to $\partial\Omega$. The scalar-valued coefficient
$\kappa$ is the permeability of the porous medium which is highly
heterogeneous in practical subsurface flow problems. By applying the
divergence theorem to Eq. 1, it is easy to validate that the source function
$f$ should satisfy the so-called compatibility condition
$\int_{\Omega}f\,\mathrm{d}\bm{x}=0$. Moreover, to guarantee the uniqueness of
the solution, an additional restriction $\int_{\Omega}p\,\mathrm{d}\bm{x}=0$
should be imposed. In reservoir simulation, the flux variable

$\bm{v}=-\kappa\nabla p,$

is usually the quantity of interest. Then we can transform the second-order
form Eq. 1 into a first-order form:

(2) $\left\\{\begin{aligned} \kappa^{-1}\bm{v}+\nabla
p&=\bm{0}\quad\text{in}\quad\Omega,\quad&&\text{(Darcy's law)},\\\
\operatorname{div}(\bm{v})&=f\quad\text{in}\quad\Omega,\quad&&\text{(mass
conservation)},\\\
\bm{v}\cdot\bm{n}&=0\quad\text{on}\quad\partial\Omega,\quad&&\text{(no-flux
boundary condition)}.\end{aligned}\right.$

We choose the no-flux boundary condition because it is much more widely
applied in reservoir simulation comparing with other types of boundary
conditions.

We define

$\bm{V}\coloneqq\left\\{\bm{v}\in
L^{2}(\Omega)^{d}\mid\operatorname{div}(\bm{v})\in L^{2}(\Omega)\ \text{and}\
\bm{v}\cdot\bm{n}=0\;\mbox{on}\;\partial\Omega\right\\},$

and

$W\coloneqq\left\\{q\in
L^{2}(\Omega)\mid\int_{\Omega}q\,\mathrm{d}\bm{x}=0\right\\}.$

Then, multiplying Darcy’s law by a test function $\bm{w}\in\bm{V}$ and the
mass conservation equation by a $q\in W$, the weak form of Eq. 2 is formulated
by seeking $(\bm{v},p)\in\bm{V}\times W$ satisfying the equations:

(3)
$\displaystyle\int_{\Omega}\kappa^{-1}\bm{v}\cdot\bm{w}\,\mathrm{d}\bm{x}-\int_{\Omega}\operatorname{div}(\bm{w})p\,\mathrm{d}\bm{x}$
$\displaystyle=0,$ $\displaystyle\forall\bm{w}\in\bm{V},$
$\displaystyle-\int_{\Omega}\operatorname{div}(\bm{v})q\,\mathrm{d}\bm{x}$
$\displaystyle=-\int_{\Omega}fq\,\mathrm{d}\bm{x},$ $\displaystyle\forall q\in
W.$

Let $\bm{V}_{h}\subset\bm{V}$ and $W_{h}\subset W$ be two finite element
spaces associated with a prescribed triangulation $\mathcal{T}_{h}$ of
$\Omega$. Then the discrete problem becomes finding
$(\bm{v}_{h},p_{h})\in\bm{V}_{h}\times W_{h}$ such that

(4)
$\displaystyle\int_{\Omega}\kappa^{-1}\bm{v}_{h}\cdot\bm{w}_{h}\,\mathrm{d}\bm{x}-\int_{\Omega}\operatorname{div}(\bm{w}_{h})p_{h}\,\mathrm{d}\bm{x}$
$\displaystyle=0,$ $\displaystyle\forall\bm{w}_{h}\in\bm{V}_{h},$
$\displaystyle-\int_{\Omega}\operatorname{div}(\bm{v}_{h})q_{h}\,\mathrm{d}\bm{x}$
$\displaystyle=-\int_{\Omega}fq_{h}\,\mathrm{d}\bm{x},$ $\displaystyle\forall
q_{h}\in W_{h}.$

There are a number of well-known families of mixed finite element spaces [7]
of $\bm{V}_{h}\times W_{h}$, here we use the lowest-order Raviart-Thomas
finite element spaces (cf. [9]) which is denoted by $\text{RT}_{0}$ in this
paper. We express the solution $(\bm{v}_{h},p_{h})$ as

$\bm{v}_{h}=\sum_{e\in\mathcal{E}^{0}_{h}}v_{e}\bm{\phi}_{e}\quad\text{and}\quad
p_{h}=\sum_{\tau\in\mathcal{T}_{h}}p_{\tau}q_{\tau},$

where $\mathcal{E}^{0}_{h}$ is the set of internal edges/faces,
$\left\\{\bm{\phi}_{e}\right\\}_{e\in\mathcal{E}^{0}_{h}}$ and
$\left\\{q_{\tau}\right\\}_{\tau\in\mathcal{T}_{h}}$ are the
bases111Generally, those bases should depend on $h$, and we drop it here for
simplicity. of $\bm{V}_{h}$ and $W_{h}$. Then algebraic representation of Eq.
4 is

(5) $\begin{pmatrix}\mathsf{M}&-\mathsf{B}^{\intercal}\\\
-\mathsf{B}&\mathsf{0}\\\ \end{pmatrix}\begin{pmatrix}\mathsf{v}\\\
\mathsf{p}\end{pmatrix}=\begin{pmatrix}0\\\ -\mathsf{f}\end{pmatrix},$

where $\mathsf{M}$ is a symmetric, positive definite matrix, $\mathsf{v}$,
$\mathsf{p}$ and $\mathsf{f}$ are column vectors.

###  2.2 The velocity elimination technique

Designing an efficient preconditioner for solving the saddle point system Eq.
5 is not an easy task. Fortunately, it is shown in [5] that if the trapezoidal
quadrature rule is applied to compute the integration
$\int_{\Omega}\kappa^{-1}\bm{v}_{h}\cdot\bm{w}_{h}\,\mathrm{d}\bm{x}$ on each
element rather than the accurate rule, then we could obtain a diagonal
velocity mass matrix $\tilde{\mathsf{M}}$ that can be inverted
straightforwardly. Therefore, the unknown $\mathsf{v}$ in Eq. 4 can be
eliminated and Eq. 5 boils down to

(6) $\mathsf{A}\mathsf{p}=\mathsf{f},$

where
$\mathsf{A}\coloneqq\mathsf{B}\tilde{\mathsf{M}}^{-1}\mathsf{B}^{\intercal}$
is a symmetric and semi-positive definite matrix. We call this technique as
velocity elimination, the equivalence between velocity elimination and a cell-
centered finite difference scheme in rectangular meshes can be referred in
[5].

Since the velocity elimination technique is closely related with the
construction of our preconditioner, we will restate it in detail here. For
simplicity, we consider a rectangular $(d=2)$ domain
$\Omega=(0,L_{x})\times(0,L_{y})$222In previous equations $\bm{x}$ is the
vector-valued spatial variable (such as Eqs. 3 and 4), while here we also use
subscripts $x$ and $y$ for $x$\- and $y$-direction. with a uniform structured
mesh $M_{x}\times M_{y}$, which says that each element is of the size of
$h_{x}\times h_{y}$ with $h_{x}=L_{x}/M_{x}$ and $h_{y}=L_{y}/M_{y}$. Note
that this mesh should be at the finest resolution and $\kappa$ should be
element-wisely constant. For an element $\tau=(i,i+1)h_{x}\times(j,j+1)h_{y}$
with $0\leq i<M_{x}$ and $0\leq j<M_{y}$, the trapezoidal quadrature rule
gives

$\displaystyle\int_{\tau}\kappa^{-1}\bm{v}_{h}\cdot\bm{w}_{h}\,\mathrm{d}\bm{x}$
$\displaystyle\approx$
$\displaystyle\frac{\left(\kappa|_{\tau}\right)^{-1}}{4}\left\\{(\bm{v}_{h}\cdot\bm{w}_{h})(p_{1})+(\bm{v}_{h}\cdot\bm{w}_{h})(p_{2})+(\bm{v}_{h}\cdot\bm{w}_{h})(p_{3})+(\bm{v}_{h}\cdot\bm{w}_{h})(p_{4})\right\\}h_{x}h_{y},$

where $\kappa|_{\tau}$ is the constant value of $\kappa$ in $\tau$,
$p_{1},\dots,p_{4}$ are the four corner points of $\tau$. According to the
definition of $\bm{\phi}_{e}$ (cf. [9]), for $p_{1}=(ih_{x},jh_{y})$, we have
$(\bm{v}_{h}\cdot\bm{w}_{h})(p_{1})=v_{e_{1}}w_{e_{1}}+v_{e_{2}}w_{e_{2}}$,
where $e_{1}$ and $e_{2}$ are the edges $\\{ih_{x}\\}\times(j,j+1)h_{y}$ and
$(i,i+1)h_{x}\times\\{jh_{y}\\}$ respectively. We can handle
$(\bm{v}_{h}\cdot\bm{w}_{h})(p_{2})$, $(\bm{v}_{h}\cdot\bm{w}_{h})(p_{3})$ and
$(\bm{v}_{h}\cdot\bm{w}_{h})(p_{4})$ similarly and then obtain an
approximation of
$\int_{\Omega}\kappa^{-1}\bm{v}_{h}\cdot\bm{w}_{h}\,\mathrm{d}\bm{x}$ as
$h_{x}h_{y}\sum_{e\in\mathcal{E}_{h}^{0}}\kappa_{e}^{-1}v_{e}w_{e}$ with
$\kappa_{e}$ is the _harmonic average_ of $\kappa_{e,+}$ and $\kappa_{e,-}$,
where $\kappa_{e,+}$ and $\kappa_{e,-}$ are two values of $\kappa$ in the two
adjacent elements along the internal edge $e$. Note that
$\int_{\Omega}\kappa^{-1}\bm{v}_{h}\cdot\bm{w}_{h}\,\mathrm{d}\bm{x}\approx
h_{x}h_{y}\sum_{e\in\mathcal{E}_{h}^{0}}\kappa_{e}^{-1}v_{e}w_{e}$ actually
provides us the expression of $\tilde{\mathsf{M}}$. For
$\int_{\Omega}\operatorname{div}(\bm{w}_{h})p_{h}\,\mathrm{d}\bm{x}$,
recalling that $p_{h}$ is also element-wisely constant and applying the
divergence theorem for each element, we have

$\displaystyle\int_{\Omega}\operatorname{div}(\bm{w}_{h})p_{h}\,\mathrm{d}\bm{x}$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}\int_{\tau}\operatorname{div}(\bm{w}_{h})p_{h}\,\mathrm{d}\bm{x}$
$\displaystyle=\sum_{\tau\in\mathcal{T}_{h}}p_{h}|_{\tau}\int_{\tau}\operatorname{div}(\bm{w}_{h})\,\mathrm{d}\bm{x}=\sum_{\tau\in\mathcal{T}_{h}}p_{h}|_{\tau}\int_{\partial\tau}\bm{w}_{h}\cdot\bm{n}\,\mathrm{d}\sigma$
$\displaystyle=-\sum_{e\in\mathcal{E}_{h}^{0}}\lvert e\rvert\llbracket
p_{h}\rrbracket_{e}w_{e},$

where $\lvert e\rvert$ is the length of $e$ and $\llbracket
p_{h}\rrbracket_{e}$ is the jump of $p_{h}$ across $e$ towards the positive
$x$\- or $y$-direction. We hence cancel out $w_{e}$ via
$h_{x}h_{y}\sum_{e\in\mathcal{E}_{h}^{0}}\kappa_{e}^{-1}v_{e}w_{e}+\sum_{e\in\mathcal{E}_{h}^{0}}\lvert
e\rvert\llbracket p_{h}\rrbracket_{e}w_{e}=0$ and obtain
$v_{e}=-\kappa_{e}\lvert e\rvert\llbracket p_{h}\rrbracket_{e}/(h_{x}h_{y})$.
Moreover, combining
$\int_{\Omega}\operatorname{div}(\bm{v}_{h})q_{h}\,\mathrm{d}\bm{x}=-\sum_{e\in\mathcal{E}_{h}^{0}}\lvert
e\rvert\llbracket q_{h}\rrbracket_{e}v_{e}$, we could derive a variational
form expression of Eq. 6 as

(7) $\sum_{e\in\mathcal{E}_{h}^{0}}\kappa_{e}\llbracket
p_{h}\rrbracket_{e}\llbracket q_{h}\rrbracket_{e}\frac{\lvert
e\rvert^{2}}{h_{x}h_{y}}=\int_{\Omega}fq_{h}\,\mathrm{d}\bm{x},\quad\forall
q_{h}\in W_{h}.$

For $d=3$, the summation $\sum_{e\in\mathcal{E}_{h}^{0}}$ in Eq. 7 should be
interpreted as the sum over all internal element faces, and
$\llbracket\cdot\rrbracket_{e}$ represents the jump across the face $e$. In
contrast to the no-flux boundary condition, for zero-Dirichlet boundary value
problems, the velocity across the boundary is not known priorly. Therefore,
following the velocity technique, additional terms will be included in Eq. 7.
Although the precise formulation is omitted here for the sake of simplicity.

Once the pressure unknowns $\mathsf{p}$ are available by solving Eq. 6, the
velocity unknowns $\mathsf{v}$ can be recovered via
$\mathsf{v}=\tilde{\mathsf{M}}^{-1}\mathsf{B}^{\intercal}\mathsf{p}$. We
emphasize that the velocity field $\bm{v}_{h}$ computed here is element-wisely
mass conservative, which is crucial for long-time subsurface flow modeling.

###  2.3 Two-level overlapping preconditioners

In this subsection, the idea of two-level overlapping preconditioners (ref.
[47]) for solving Eq. 6 is briefly presented. Assume the domain $\Omega$ is
divided into a union of _disjoint_ polygonal subdomains denoted by
$\left\\{K_{i}\right\\}_{i=1}^{N}$ with each $K_{i}$ consisting of fine
elements from $\mathcal{T}_{h}$. Therefore, we obtain a coarse “triangulation”
of $\Omega$ that is denoted by $\mathcal{T}_{H}$, and we will also call a
subdomain $K_{i}$ as a coarse element. Let
$\left\\{K_{i}^{m}\right\\}_{i=1}^{N}$ be an overlapping partition by
extending $m$ layers of fine elements to each coarse element $K_{i}$. In later
presentation, we also name $K_{i}^{m}$ as an oversampling coarse element, and
the term “oversampling” is coined in the oversampling multiscale finite
element method [33]. Figure 1 is an illustration of the two-scale mesh, a fine
element $\tau$, a coarse element $K_{i}$ and its oversampling coarse element
$K_{i}^{m}$ with $m=2$.

$K_{i}$$K_{i}^{2}$$\tau$ Figure 1: An illustration of the two-scale mesh, a
fine element $\tau$, a coarse element $K_{i}$ and its oversampling coarse
element $K_{i}^{m}$ with $m=2$.

We introduce coarse basis functions
$\left\\{\Phi_{i}\right\\}_{i=1}^{N^{\text{c}}}$ associated with the coarse
mesh $\mathcal{T}_{H}$, where $N^{\text{c}}$ is the number of bases, and
$N^{\text{c}}$ in generally is different from the number of coarse elements
$N$. Note that each basis should be resolved on the fine mesh
$\mathcal{T}_{h}$, and we can hence represent the coarse space as

$W_{H}^{\text{c}}\coloneqq\operatorname{span}\left\\{\Phi_{i}\right\\}_{i=1}^{N^{\text{c}}}$

and its projection to the fine space as a prolongation matrix
$\mathsf{R}_{0}^{\intercal}$ with each column of $\mathsf{R}_{0}^{\intercal}$
corresponding to a coarse basis. The coarse space can simply be
$\text{RT}_{0}$ space for pressure on the coarse mesh $\mathcal{T}_{H}$, it
can also be a _generalized multiscale finite element space_. After determining
$W_{H}^{\text{c}}$, the matrix of the global solver can be obtained by
$\mathsf{A}_{0}\coloneqq\mathsf{R}_{0}\mathsf{A}\mathsf{R}_{0}^{\intercal}$,
where the size of $\mathsf{A}_{0}$ is $N^{\text{c}}\times N^{\text{c}}$. We
emphasize that constructing a low-dimensional but effective coarse space is
pivotal in designing two-level preconditioners.

An indispensable component in two-level preconditioners is so-called local
solvers. In our paper, the local solvers are built from discretizing zero-
Dirichlet problems defined on _oversampling_ coarse elements
$\left\\{K_{i}^{m}\right\\}_{i=1}^{N}$. That is, we solve the following
equation on the fine mesh $\mathcal{T}_{h}$ via the $\text{RT}_{0}$ method:

(8) $\left\\{\begin{aligned} \kappa^{-1}\bm{v}+\nabla
p&=\bm{0}\quad&&\text{in}\quad K_{i}^{m},\\\
\nabla\cdot\bm{v}&=f\quad&&\text{in}\quad K_{i}^{m},\\\
p&=0\quad&&\text{on}\quad\partial K_{i}^{m}\cap\Omega,\\\
\bm{v}\cdot\bm{n}&=0\quad&&\text{on}\quad\partial
K_{i}^{m}\cap\partial\Omega.\end{aligned}\right.$

Note that we utilize the velocity elimination technique again to obtain a
linear system that is only related to pressure. We denote $\mathsf{A}_{i}$ the
matrix of the local solver corresponding to $K_{i}^{m}$, and $\mathsf{R}_{i}$
the restriction matrix that restricts the Degrees of Freedom (DoF) of pressure
from $\Omega$ to $K_{i}^{m}$. Note that $\mathsf{A}_{i}$ is not equal to
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$, and this fact is
different from the discretization by the Lagrange finite element method.
Except for special situations (e.g., $\partial
K_{i}^{m}\cap\Omega=\varnothing$), the matrix $\mathsf{A}_{i}$ is invertible,
and we will always assume that $\mathsf{A}_{i}^{-1}$ exists.

The fine-scale linear system Eq. 6 will be solved by preconditioned iterative
solvers such as Generalized Minimal Residual (GMRES) methods with a two-level
overlapping preconditioner of the form:

$\mathsf{P}^{-1}=\mathsf{R}_{0}^{\intercal}\mathsf{A}_{0}^{\dagger}\mathsf{R}_{0}+\sum_{i=1}^{N}\mathsf{R}_{i}^{\intercal}\mathsf{A}_{i}^{-1}\mathsf{R}_{i},$

where $\mathsf{A}_{0}^{\dagger}$ is the pseudoinverse of $\mathsf{A}_{0}$
because $\mathsf{A}_{0}$ is generally singular and “$-1$” of $\mathsf{P}^{-1}$
is just a conventional notation. The definition of this preconditioner
indicates that it needs to solve a global system and multiple local problems
in each iteration. In practice, we will factorize $\mathsf{A}_{0}$ and
$\mathsf{A}_{i}$ during the preparation phase, which saves computing time on
solving precondition systems in the latter iteration phase.

##  3 GMsFEM based two-level overlapping preconditioner

The choice of the coarse space $W_{H}^{\text{c}}$ has a significant influence
on the performance, and we will adopt the methodology of GMsFEMs to construct
a coarse space that incorporates heterogeneity information of $\kappa$. Follow
the simplification in Section 2.2, we will still consider a rectangular domain
with a uniform structured mesh, and the generalization to unstructured meshes
is also possible (see [28]). Meanwhile, the coarse mesh $\mathcal{T}_{H}$ is
also structured and each coarse element $K_{i}$ is a rectangular subdomain.

In GMsFEMs [21], two bilinear forms $a_{i}(\cdot,\cdot)$ and
$s_{i}(\cdot,\cdot)$ corresponding to the $i$-th subdomain333Subdomains
$\left\\{\omega_{i}\right\\}$ are overlapped in [21], while coarse elements in
our paper are all disjointed. are important in selecting effective low-
dimensional representations from snapshot spaces. According to Eq. 7, we
define $a_{i}(\cdot,\cdot)$ as

(9)
$a_{i}(q_{h}^{\prime},q_{h})\coloneqq\sum_{e\in\mathcal{E}_{h}^{0}(K_{i})}\kappa_{e}\llbracket
q_{h}^{\prime}\rrbracket_{e}\llbracket q_{h}\rrbracket_{e}\frac{\lvert
e\rvert^{2}}{h_{x}h_{y}},\quad\forall q_{h},q_{h}^{\prime}\in W_{h}(K_{i}),$

where $\mathcal{E}_{h}^{0}(K_{i})$ is the internal edge set in the coarse
element $K_{i}$ and $W_{h}(K_{i})$ is the restriction space of $W_{h}$ on
$K_{i}$. Due to that our discretization is not a Galerkin method, the
expression of $a_{i}(\cdot,\cdot)$ is different from the original one in [21].
We define $s_{i}(\cdot,\cdot)$ as

(10)
$s_{i}(q_{h}^{\prime},q_{h})\coloneqq\int_{K_{i}}\tilde{\kappa}q_{h}^{\prime}q_{h}\,\mathrm{d}\bm{x},\quad\forall
q_{h},q_{h}^{\prime}\in W_{h}(K_{i}),$

where $\tilde{\kappa}$ must satisfy Assumption B as described in Section 4. In
addition, we explain in Section 5 that it is also reasonable to set
$\tilde{\kappa}=\kappa$. Then, we solve the following spectral problem in each
coarse element $K_{i}$:

(11)
$\sum_{e\in\mathcal{E}_{h}^{0}(K_{i})}\kappa_{e}\llbracket\Phi_{h}\rrbracket_{e}\llbracket
q_{h}\rrbracket_{e}\frac{\lvert
e\rvert^{2}}{h_{x}h_{y}}=\lambda\int_{K_{i}}\tilde{\kappa}\Phi_{h}q_{h}\,\mathrm{d}\bm{x},\quad\forall
q_{h}\in W_{h}(K_{i}),$

where $\Phi_{h}$ is an eigenvector corresponding to the eigenvalue $\lambda$.
Note that Eq. 11 is the discretization of the spectral problem—find $\Phi\in
H^{1}(K_{i})$ such that
$-\operatorname{div}\left(\kappa\nabla\Phi\right)=\lambda\tilde{\kappa}\Phi$
in $K_{i}$—via $\text{RT}_{0}$ elements and the velocity elimination
technique. These eigenvectors supported in coarse elements are solutions of
local spectral problems that contain fine-scale information $\kappa$,
therefore we also call them multiscale basis functions. After solving the
eigenvalue problem Eq. 11, eigenvectors $\\{\Phi_{h,i}^{j}\\}_{j=1}^{L_{i}}$
corresponding to the $L_{i}$ smallest eigenvalues will be collected to form
the local coarse space $W_{H}^{\text{c}}(K_{i})$ as
$W_{H}^{\text{c}}(K_{i})=\operatorname{span}\\{\Phi_{h,i}^{j}:j=1,\dots,L_{i}\\}$.
Then, the global coarse space $W_{H}^{c}$ is the direct sum444Two different
local coarse space are associated to two different subdomains, which says the
direct sum operation is inappropriate. Strictly, the direct sum here should
implicitly come after performing zero-extensions of bases to the whole domain.
of $W_{H}^{c}(K_{i})$ as
$W_{H}^{\text{c}}=\bigoplus_{i=1}^{N}W_{H}^{\text{c}}(K_{i})$. The dimension
of $W_{H}^{c}$ equals the total number of eigenvectors, that is
$N^{\text{c}}=\sum_{i=1}^{N}L_{i}$. We emphasize again that Eq. 11 is solved
in each coarse element $K_{i}$ that has no overlap with its neighbors, and
therefore no communication is required when adopting parallel computing.
Comparing with previous efforts [23, 24, 43, 44], in which local problems are
solved in overlapping subdomains, the advantage of our method in parallel
computing implementation is remarkable.

######  Remark 3.1.

In the original paper on GMsFEMs [21], the authors propose the use of
$\kappa$-harmonic spaces as snapshot spaces, which involves solving the
spectral problem in Eq. 11 on a subspace of $W_{h}(K_{i})$. In our paper, we
advocate for the approach in Eq. 11 for eigenvalue problems, which could be
understood by utilizing full $W_{h}(K_{i})$ as snapshot spaces. While using
$\kappa$-harmonic spaces can reduce the degrees of freedom in spectral
problems, it can also result in dense eigenvalue systems. Furthermore,
implementing the use of $\kappa$-harmonic spaces is more complex, making it
challenging to assess the actual performance improvement. Additionally, our
numerical findings indicate that in parallel computing scenarios, the time
spent on spectral problems does not significantly contribute to the overall
computational time.

The bilinear form $s_{i}(\cdot,\cdot)$ is symmetric and positive definite,
while $a_{i}(\cdot,\cdot)$ is symmetric but always semi-definite, which
implies the first eigenvalue of Eq. 11 is always $0$. Moreover, it is easy to
see there is a unique and normalized eigenvector corresponding to the
eigenvalue $0$, and this eigenvector is exactly the constant vector. Note that
the $\text{RT}_{0}$ spaces $\mathbf{V}_{H}\times W_{H}$ on the coarse mesh
$\mathcal{T}_{H}$, bases of the pressure part $W_{H}$ are also piece-wisely
constant. Then, our coarse space can be viewed as an enrichment of $W_{H}$,
which is a standard coarse space in two-grid preconditioners [10]. Recall that
the no-flux boundary condition is imposed for our model problem Eq. 1, and the
matrix $\mathsf{A}$ of the fine-scale linear system Eq. 6 is actually singular
with constant vectors as its kernel. Note that the coarse space
$W_{H}^{\text{c}}$ also contains the kernel of $\mathsf{A}$, and this implies
that $W_{H}^{\text{c}}$ inherits some physical information of the original
system and can be interpreted as an appropriate upscaling of the fine-scale
space.

Figure 2 demonstrates eigenvalues computed by setting different permeability
profiles. In the left column, two types of media (3- and 5-channel
configurations) that contain long channels are visualized, and we will fix
$\kappa\equiv 1$ in the background region (purple) and
$\kappa\equiv\kappa^{*}$ in channels (yellow), where $\kappa^{*}$ is from
$\left\\{10^{0},10^{1},10^{2},10^{3},10^{4}\right\\}$ to see the effect of
contrast ratios. We plot the $11$ smallest eigenvalues of Eq. 11 with
$\tilde{\kappa}=\kappa$ in the middle column, with $\tilde{\kappa}=1$ in the
right column, where plots in the first row are corresponding to the 3-channel
configuration and the second row is for the 5-channel configuration
accordingly. We can observe the patterns of eigenvalues for
$\tilde{\kappa}=\kappa$ and $1$ are significantly different. Checking vertical
axes of the plots in the right and middle columns, eigenvalues corresponding
to $\tilde{\kappa}=1$ seem to be sensitive to different configurations, while
eigenvalues to $\tilde{\kappa}=\kappa$ are more stable. An interesting
phenomenon could be noticed in the middle column is for the 3-channel
(5-channel) configuration, the $3$ ($5$) smallest eigenvalues will approach to
$0$ as increasing contrast ratios. Actually, it is shown mathematically that
in [23], if there are $m$ channels, the $m$ smallest eigenvalues may approach
to $0$ as increasing contrast ratios, which is not favored in building robust
preconditions. Generally, We need to include enough eigenvectors into our
coarse space, and the number of channels may provide a guidance. However, for
a complex medium, it is almost impossible to have a priori knowledge about the
precise number of channels. In practice, we may preset a number or choose
eigenvectors adaptively such that the largest eigenvalue crosses a certain
threshold.

Figure 2: (left column) the 3- and 5-channel configurations; (middle column)
the $11$ smallest eigenvalues calculated by setting $\tilde{\kappa}=\kappa$
w.r.t. different $\kappa^{*}$ that is the value of $\kappa$ in channels, and
the top (bottom) plot is corresponding to the 3-channel (5-channel)
configuration; (right column) the $11$ smallest eigenvalues calculated by
setting $\tilde{\kappa}=1$ w.r.t. different $\kappa^{*}$, and the top (bottom)
plot is corresponding to the 3-channel (5-channel) configuration.

We summarize our method of constructing the preconditioner in Algorithm 1.
Basically, the method consists of two phases—preparation and iteration. Note
that every step listed in Algorithm 1 could be accelerated in a parallel
framework. Although the exact iteration phase depends on the iterative solver,
we just focus on applying $\mathsf{P}^{-1}\mathsf{r}$, where $\mathsf{r}$ is
usually a residual vector.

Algorithm 1 GMsFEM based two-level overlapping preconditioner

1:Preparation phase:

2: Set up the decomposition $\mathcal{T}_{H}$ for the domain $\Omega$; load
the coefficient profile $\kappa$; set the source term $f$.

3: Construct the matrix $\mathsf{A}$ and the right-hand vector $\mathsf{f}$ of
Eq. 6.

4: For each oversampling coarse element $K_{i}^{m}$, construct the matrix
$\mathsf{A}_{i}$ by discretizing Eq. 8 via velocity elimination; factorize
$\mathsf{A}_{i}$.

5: For each coarse element $K_{i}$, solve the spectral problem Eq. 11 and
obtain $L_{i}$ eigenvectors; form the coarse space $W_{H}^{\text{c}}$ or the
prolong matrix $\mathsf{R}_{0}^{\intercal}$.

6: Construct the matrix $\mathsf{A}_{0}$ via the coarse space
$W_{H}^{\text{c}}$; factorize $\mathsf{A}_{0}$.

7:Iteration phase (apply $\mathsf{P}^{-1}\mathsf{r}$):

8: Restrict $\mathsf{r}$ to the oversampling coarse element $K_{i}^{m}$ and
obtain $\mathsf{r}_{i}$; use factorized $\mathsf{A}_{i}$ to get
$\mathsf{A}_{i}^{-1}\mathsf{r}_{i}$; sum all
$\mathsf{A}_{i}^{-1}\mathsf{r}_{i}$ according to DoF in the fine mesh.

9: Project $\mathsf{r}$ to the coarse space $W_{H}^{\textsc{c}}$ and obtain
$\mathsf{r}_{0}=\mathsf{R}_{0}\mathsf{r}$; use factorized $\mathsf{A}_{0}$ to
get $\mathsf{A}_{0}^{\dagger}\mathsf{r}_{0}$; project the solution
$\mathsf{A}_{0}^{\dagger}\mathsf{r}_{0}$ back to the fine space with the
coarse bases.

10: Take a summation of the resulting vectors in the aforementioned two steps
and obtain $\mathsf{P}^{-1}\mathsf{r}$.

##  4 Analysis

For iterative solvers, a smaller iteration number usually means a better
performance555Actually, the performance or the elapsed time for a
preconditioned iterative solver cannot be solely determined by numbers of
iterations. To see this, just assume the precondition system is the original
system, then only one iteration is needed while there is no acceleration
benefiting from preconditioning at all., and numbers of iterations are
dramatically affected by the condition number of the preconditioned matrix
[45]. However, the canonical definition of condition numbers is for invertible
matrices, while the matrix $\mathsf{A}$ in Eq. 6 is singular. According to
[11], for a linear system $\mathsf{M}\mathsf{x}=\mathsf{b}$, if the right-hand
vector $\mathsf{b}$ belongs to $\operatorname{im}(\mathsf{M})$ the image of
$\mathsf{M}$, then GMRES can safely find a least-squares solution or the
pseudoinverse solution. Moreover, the paper [11] suggests a redefinition of
the condition number for $\mathsf{A}$ as

$\operatorname{cond}(\mathsf{A})=\frac{\sigma_{\max}(\mathsf{A})}{\sigma_{2}(\mathsf{A})},$

where $\sigma_{\max}(\mathsf{A})$ is the largest singular value of
$\mathsf{A}$ and $\sigma_{2}(\mathsf{A})$ is the second least singular value.
Hence, this section aims to show
$\operatorname{cond}(\mathsf{P}^{-1}\mathsf{A})$ can be bounded above with
suitable assumptions.

The main tool is the so-called fictitious space lemma (ref. [26]). Before
presenting this lemma, we introduce several notations. Let $n$ be the number
of total DoF of pressure, $n_{i}$ be the number of DoF of pressure restricted
on the coarse element $K_{i}$, and $n_{i}^{m}$ be the number of DoF of
pressure restricted on the oversampling coarse element $K_{i}^{m}$,
$\underline{V}$ be the product space
$\mathbb{R}^{N^{\text{c}}}\times\mathbb{R}^{n_{1}^{m}}\times\dots\times\mathbb{R}^{n_{N}^{m}}$.
Then we have $\mathsf{A}\in\mathbb{R}^{n\times n}$, and
$\mathsf{A}_{i}\in\mathbb{R}^{n_{i}^{m}\times
n_{i}^{m}},\mathsf{R}_{i}\in\mathbb{R}^{n_{i}^{m}\times n}$ for $i=1,\dots,N$.
We define $\mathsf{a}(\cdot,\cdot)$ as a symmetric, positive semi-definite
bilinear form on $\mathbb{R}^{n}$ by

$\mathsf{a}(\mathsf{v},\mathsf{w})=\mathsf{v}\cdot\mathsf{A}\mathsf{w},\quad\forall\mathsf{v},\mathsf{w}\in\mathbb{R}^{n},$

which is exactly the left-hand section of Eq. 7 in a linear algebraic form.
Recall that the dimension of the coarse space $W_{H}^{\text{c}}$ is
$N^{\text{c}}$, which gives $\mathsf{A}_{0}\in\mathbb{R}^{N^{\text{c}}\times
N^{\text{c}}}$ and $\mathsf{R}_{0}\in\mathbb{R}^{N^{\text{c}}\times n}$. The
fictitious space lemma is stated below, and we avoid using abstract Hilbert
spaces here.

######  Lemma 4.1.

Let $\mathsf{b}(\cdot,\cdot)$ be a symmetric, positive semi-definite bilinear
form on $\underline{V}$ with a matrix representation as

$\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{w}})=\underline{\mathsf{v}}\cdot\underline{\mathsf{B}}\,\underline{\mathsf{w}},\quad\forall\underline{\mathsf{v}},\underline{\mathsf{w}}\in\underline{V},$

$\mathcal{R}:\underline{V}\rightarrow\mathbb{R}^{n}$ be a linear map, and
$\mathcal{R}^{*}:\mathbb{R}^{n}\rightarrow\underline{V}$ be an adjoint
operator by

$\mathcal{R}\underline{\mathsf{v}}\cdot\mathsf{w}=\underline{\mathsf{v}}\cdot\mathcal{R}^{*}\mathsf{w},\quad\forall\underline{\mathsf{v}}\in\underline{V},\mathsf{w}\in\mathbb{R}^{n}.$

Assume that $\mathcal{R}$ maps $\ker(\underline{\mathsf{B}})$ into
$\ker(\mathsf{A})$, and there are positive constants $C_{T}$ and $C_{R}$, such
that for all $\mathsf{u}\in\mathbb{R}^{n}$ there exists a
$\underline{\mathsf{v}}\in\underline{V}$ with
$\mathsf{u}=\mathcal{R}\underline{\mathsf{v}}$ and

(12)
$C_{T}\,\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}})\leq\mathsf{a}(\mathsf{u},\mathsf{u}),$

and also

(13)
$\mathsf{a}(\mathcal{R}\underline{\mathsf{w}},\mathcal{R}\underline{\mathsf{w}})\leq
C_{R}\,\mathsf{b}(\underline{\mathsf{w}},\underline{\mathsf{w}}),\quad\forall\underline{\mathsf{w}}\in\underline{V}.$

Then for all $\mathsf{u}\in\mathbb{R}^{n}$ with
$\mathsf{a}(\mathsf{u},\mathsf{u})>0$,

$C_{T}\,\mathsf{a}(\mathsf{u},\mathsf{u})\leq\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})\leq
C_{R}\,\mathsf{a}(\mathsf{u},\mathsf{u}),$

where $\underline{\mathsf{B}}^{\dagger}$ is the pseudoinverse of
$\underline{\mathsf{B}}$.

The proof of Lemma 4.1 is an extension of the proof in [26], while several
modifications required due to that $\mathsf{a}(\cdot,\cdot)$ and
$\mathsf{b}(\cdot,\cdot)$ are no longer positive definite.

######  Proof 4.2.

Note that
$\underline{\mathsf{v}}-\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{v}}\in\ker(\underline{\mathsf{B}})$
for all $\underline{\mathsf{v}}\in\underline{V}$. Thanks to the assumption
that $\mathcal{R}$ maps $\ker(\underline{\mathsf{B}})$ into
$\ker(\mathsf{A})$, we can show that

$\mathcal{R}\underline{\mathsf{v}}\cdot\mathsf{A}\mathsf{w}=\mathcal{R}\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{v}}\cdot\mathsf{A}\mathsf{w}$

for all $\mathsf{w}\in\mathbb{R}^{n}$. Then, for all
$\mathsf{u}\in\mathbb{R}^{n}$ with $\mathsf{a}(\mathsf{u},\mathsf{u})>0$,
choosing $\underline{\mathsf{v}}\in\underline{V}$ such that
$\mathsf{u}=\mathcal{R}\underline{\mathsf{v}}$, we have

$\displaystyle\mathsf{a}(\mathsf{u},\mathsf{u})=$
$\displaystyle\mathsf{u}\cdot\mathsf{A}\mathsf{u}=\mathcal{R}\underline{\mathsf{v}}\cdot\mathsf{A}\mathsf{u}=\mathcal{R}\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{v}}\cdot\mathsf{A}\mathsf{w}=\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{v}}\cdot\mathcal{R}^{*}\mathsf{A}\mathsf{u}=\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{v}}\cdot\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}$
$\displaystyle\leq$
$\displaystyle\left(\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{v}}\cdot\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{v}}\right)^{1/2}\left(\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\sqrt{\underline{\mathsf{B}}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\right)^{1/2}$
$\displaystyle=$
$\displaystyle\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}})^{1/2}\left(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\mathsf{A}\mathsf{u}\right)^{1/2}$
$\displaystyle\leq$ $\displaystyle
C_{T}^{-1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}\left(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\mathsf{A}\mathsf{u}\right)^{1/2}$
$\displaystyle=$ $\displaystyle
C_{T}^{-1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})^{1/2},$

which gives that
$C_{T}\,\mathsf{a}(\mathsf{u},\mathsf{u})\leq\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})$
the first part of the target inequality. For the second part,

$\displaystyle\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})=$
$\displaystyle\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\mathsf{A}\mathsf{u}=\sqrt{\mathsf{A}}\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\sqrt{\mathsf{A}}\mathsf{u}$
$\displaystyle\leq$
$\displaystyle\left(\sqrt{\mathsf{A}}\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\sqrt{\mathsf{A}}\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\right)^{1/2}\left(\sqrt{\mathsf{A}}\mathsf{u}\cdot\sqrt{\mathsf{A}}\mathsf{u}\right)^{1/2}$
$\displaystyle=$
$\displaystyle\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u})^{1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}$
$\displaystyle\leq$ $\displaystyle
C_{R}^{1/2}\mathsf{b}(\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u})^{1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}$
$\displaystyle=$ $\displaystyle
C_{R}^{1/2}\left(\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\underline{\mathsf{B}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\right)^{1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}$
$\displaystyle=$ $\displaystyle
C_{R}^{1/2}\left(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\underline{\mathsf{B}}\,\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u}\cdot\mathsf{A}\mathsf{u}\right)^{1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2}$
$\displaystyle=$ $\displaystyle
C_{R}^{1/2}\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})^{1/2}\mathsf{a}(\mathsf{u},\mathsf{u})^{1/2},$

and this gives
$\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})\leq
C_{R}\,\mathsf{a}(\mathsf{u},\mathsf{u})$ due to
$0<C_{T}\,\mathsf{a}(\mathsf{u},\mathsf{u})\leq\mathsf{a}(\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{*}\mathsf{A}\mathsf{u},\mathsf{u})$.

If take $\mathsf{b}(\cdot,\cdot)$ as

(14)
$\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{w}})=\underline{\mathsf{v}}\cdot\underline{\mathsf{B}}\,\underline{\mathsf{w}}\coloneqq\mathsf{v}_{0}\cdot\mathsf{A}_{0}\mathsf{w}_{0}+\sum_{i=1}^{N}\mathsf{v}_{i}\cdot\mathsf{A}_{i}\mathsf{w}_{i},$

where
$\underline{\mathsf{v}}=[\mathsf{v}_{0},\mathsf{v}_{1},\dots,\mathsf{v}_{N}],\underline{\mathsf{w}}=[\mathsf{w}_{0},\mathsf{w}_{1},\dots,\mathsf{w}_{N}]\in\underline{V}$,
and the map $\mathcal{R}$ as

(15)
$\mathcal{R}\underline{\mathsf{v}}\coloneqq\mathsf{R}_{0}^{\intercal}\mathsf{v}_{0}+\sum_{i=1}^{N}\mathsf{R}_{i}^{\intercal}\mathsf{v}_{i},$

we can then show that

$\underline{\mathsf{B}}^{\dagger}\,\underline{\mathsf{v}}=[\mathsf{A}_{0}^{\dagger}\mathsf{v}_{0},\mathsf{A}_{1}^{-1}\mathsf{v}_{1},\dots,\mathsf{A}_{N}^{-1}\mathsf{v}_{N}],\quad\forall\underline{\mathsf{v}}\in\underline{V},$

and

$\mathcal{R}^{*}\mathsf{v}=[\mathsf{R}_{0}\mathsf{v},\mathsf{R}_{1}\mathsf{v},\dots,\mathsf{R}_{N}\mathsf{v}],\quad\forall\mathsf{v}\in\mathbb{R}^{n}.$

It is obvious to examine that
$\mathcal{R}\underline{\mathsf{B}}^{\dagger}\mathcal{R}^{\star}\mathsf{A}=\mathsf{P}^{-1}\mathsf{A}$,
and Lemma 4.1 tells

$\operatorname{cond}(\mathsf{P}^{-1}\mathsf{A})\leq\frac{C_{T}}{C_{R}}$

given that Eqs. 12 and 13 hold, which boils our objective down to estimations
of $C_{T}$ and $C_{R}$.

We have mentioned that
$\mathsf{A}_{i}\neq\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$.
Actually, the relation could be further elucidated as the following lemma,
where the notation $\mathsf{M}\lesssim\mathsf{N}$ for two matrices
$\mathsf{M}$ and $\mathsf{N}$ means $\mathsf{N}-\mathsf{M}$ is positive semi-
definite.

######  Lemma 4.3.

The relation
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}\lesssim\mathsf{A}_{i}$
holds.

$\mathtt{k}_{p,q}$$\mathtt{k}_{p-1,q}$$\mathtt{k}_{p+1,q}$$\mathtt{k}_{p,q-1}$$\mathtt{k}_{p,q+1}$$K_{i}^{2}$$K_{i}$
Figure 3: An illustration for the proof of Lemma 4.3, where
$\mathtt{k}_{p,q}$, $\mathtt{k}_{p-1,q}$, $\mathtt{k}_{p+1,q}$,
$\mathtt{k}_{p,q-1}$ and $\mathtt{k}_{p,q+1}$ are values of $\kappa$ on the
fine elements respectively, the oversampling coarse element $K_{i}^{m}$ with
$m=2$ is filled with gray and also decorated with dashed borderlines, and the
fine element in the top-left corner of $K_{i}^{m}$ is highlighted.

######  Proof 4.4.

For a fine element $\tau\in\mathcal{T}_{h}$ with $\tau\subset K_{i}^{m}$,
which is corresponding to the $s$-th entry of a vector in
$\mathbb{R}^{n_{i}^{m}}$, the $s$-th row of
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ is exactly same with the
$s$-th row of $\mathsf{A}_{i}$ if $\tau$ does not contact with $\partial
K_{i}^{m}$. This is because $\mathsf{A}_{i}$ is from discretizing Eq. 8 which
shares the same law and discretization of the original system Eq. 2 except for
the boundary condition. Therefore, we only need to check entries of
$\mathsf{A}_{i}$ and $\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ that
are related to fine elements which contact the boundary.

The derivation will be easier to present if we consider a $2$D uniform
structured mesh as in Section 2.2. Moreover, we take the fine element in the
top-left corner of $K_{i}^{m}$ to study the corresponding rows of
$\mathsf{A}_{i}$ and $\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$. We
illustrate the geometric information in Fig. 3, where $\mathtt{k}_{p,q}$,
$\mathtt{k}_{p-1,q}$, $\mathtt{k}_{p+1,q}$, $\mathtt{k}_{p,q-1}$ and
$\mathtt{k}_{p,q+1}$ are values of $\kappa$ on the fine elements respectively.
Let $s$ be the index of the entry related to the top-left fine element, and we
can see that the $s$-th row of $\mathsf{A}_{i}$ differs from the $s$-th row of
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ only at
$[\mathsf{A}_{i}]_{s,s}$ and
$[\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}]_{s,s}$, which can be
calculated as

$\displaystyle{}_{s,s}=$
$\displaystyle\left(\frac{2}{1/\mathtt{k}_{p,q}}+\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p+1,q}}\right)\frac{h_{y}}{h_{x}}$
$\displaystyle+\left(\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p,q-1}}+\frac{2}{1/\mathtt{k}_{p,q}}\right)\frac{h_{x}}{h_{y}},$
$\displaystyle[\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}]_{s,s}=$
$\displaystyle\left(\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p-1,q}}+\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p+1,q}}\right)\frac{h_{y}}{h_{x}}$
$\displaystyle+\left(\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p,q-1}}+\frac{2}{1/\mathtt{k}_{p,q}+1/\mathtt{k}_{p,q+1}}\right)\frac{h_{x}}{h_{y}}.$

We can hence see that
$[\mathsf{A}_{i}]_{s,s}\leq[\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}]_{s,s}$.
From this case, we can see that
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}-\mathsf{A}_{i}$ is exactly
a diagonal matrix that all entries are non-negative, which essentially says
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}\lesssim\mathsf{A}_{i}$.

######  Remark 4.5.

The proof of Lemma 4.3 also explains why we use $\mathsf{A}_{i}$ rather that
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ for local solvers.
Because constructing $\mathsf{A}_{i}$ needs the values of $\kappa$ in
$K_{i}^{m}$, while constructing
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ needs the values of
$\kappa$ in $K_{i}^{m+1}$. In a message-passing model of parallel computing,
each process handles a connected subdomain that consists of several coarse
elements from $\\{K_{i}\\}_{i=1}^{N}$ and $m$ layers of ghost points for
inter-communications. For
$\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}$, we need to set up $m+1$
layers, which means more communications and will largely degrade the overall
performance.

We need the first assumption about the oversampling partition
$\\{K_{i}^{m}\\}_{i=1}^{N}$.

#### Assumption A

Suppose that $m\geq 1$ and there exists a positive constant $C_{\text{os}}$
such that for any $i$

$\operatorname{card}\left\\{K_{j}^{m}:K_{j}^{m}\cap
K_{i}^{m+1}\neq\varnothing,1\leq j\leq N\right\\}\leq C_{\text{os}}.$

Note that $\mathsf{R}_{j}\mathsf{A}\mathsf{R}_{i}^{\intercal}$ is a zero
matrix if $K_{j}^{m}\cap K_{i}^{m+1}=\varnothing$ and $1\leq i,j\leq N$, while
$K_{j}^{m}\cap K_{i}^{m}=\varnothing$ cannot result in
$\mathsf{R}_{j}\mathsf{A}\mathsf{R}_{i}^{\intercal}=\mathsf{0}$. The next
lemma provides an estimation of $C_{R}$ in Eq. 13.

######  Lemma 4.6.

Under assumption A, for all $\underline{\mathsf{v}}\in\underline{V}$,

$\mathsf{a}(\mathcal{R}\underline{\mathsf{v}},\mathcal{R}\underline{\mathsf{v}})\leq\left(1+4C_{\text{os}}\right)\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}}),$

where $\mathcal{R}$ and $\mathsf{b}(\cdot,\cdot)$ are defined in Eqs. 15 and
14 respectively.

######  Proof 4.7.

The inequality $\mathsf{a}(\mathsf{v}+\mathsf{w},\mathsf{v}+\mathsf{w})\leq
2\left(\mathsf{a}(\mathsf{v},\mathsf{v})+\mathsf{a}(\mathsf{w},\mathsf{w})\right)$
always holds if $\mathsf{A}$ is positive semi-definite. Then we have,

$\mathsf{a}(\mathcal{R}\underline{\mathsf{v}},\mathcal{R}\underline{\mathsf{v}})\leq
2\left\\{\mathsf{a}(\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0},\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0})+\mathsf{a}(\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i},\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i})\right\\}.$

Let $J_{i}\coloneqq\\{j:K_{j}^{m}\cap K_{i}^{m+1}\neq\varnothing,1\leq j\neq
N\\}$. We can see that $\operatorname{card}(J_{i})\leq C_{\text{os}}$
according to Assumption A and

$\displaystyle\mathsf{a}(\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i},\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i})=$
$\displaystyle\sum_{i,j=1}^{N}\mathsf{v}_{i}\cdot\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{j}^{\intercal}\mathsf{v}_{j}=\sum_{i=1}^{N}\sum_{j\in
J_{i}}\mathsf{v}_{i}\cdot\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{j}^{\intercal}\mathsf{v}_{j}$
$\displaystyle\leq$ $\displaystyle 2\sum_{i=1}^{N}\sum_{j\in
J_{i}}\mathsf{v}_{i}\cdot\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}\mathsf{v}_{i}+\mathsf{v}_{j}\cdot\mathsf{R}_{j}\mathsf{A}\mathsf{R}_{j}^{\intercal}\mathsf{v}_{j}$
$\displaystyle=$ $\displaystyle
2C_{\text{os}}\sum_{i=1}^{N}\mathsf{v}_{i}\cdot\mathsf{R}_{i}\mathsf{A}\mathsf{R}_{i}^{\intercal}\mathsf{v}_{i}+2\sum_{i=1}^{N}\sum_{j\in
J_{i}}\mathsf{v}_{j}\cdot\mathsf{R}_{j}\mathsf{A}\mathsf{R}_{j}^{\intercal}\mathsf{v}_{j}$
$\displaystyle\leq$ $\displaystyle
2C_{\text{os}}\sum_{i=1}^{N}\mathsf{v}_{i}\cdot\mathsf{A}_{i}\mathsf{v}_{i}+2\sum_{i=1}^{N}\sum_{j\in
J_{i}}\mathsf{v}_{j}\cdot\mathsf{A}_{j}\mathsf{v}_{j},$

where the last line is from Lemma 4.3. Note that in the summation
$\sum_{i=1}^{N}\sum_{j\in J_{i}}$, each index $i\in\\{1,\dots,N\\}$ will be
counted at most $C_{\text{os}}$ times. Therefore, we obtain

$\mathsf{a}(\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i},\sum_{i=1}^{N}\mathsf{R}^{\intercal}_{i}\mathsf{v}_{i})\leq
4C_{\text{os}}\sum_{i=1}^{N}\mathsf{v}_{i}\cdot\mathsf{A}_{i}\mathsf{v}_{i}.$

Recalling that
$\mathsf{a}(\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0},\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0})=\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0}\cdot\mathsf{A}\mathsf{R}^{\intercal}_{0}\mathsf{v}_{0}=\mathsf{v}_{0}\cdot\mathsf{A}_{0}\mathsf{v}_{0}$,
we hence finish the proof.

We can represent $s_{i}(\cdot,\cdot)$ defined in Eq. 10 in a linear algebraic
form as

$\mathsf{s}_{i}(\mathsf{v},\mathsf{w})=\mathsf{v}\cdot\mathsf{S}_{i}\mathsf{w}$

for all $\mathsf{v}$ and $\mathsf{w}$ in $\mathbb{R}^{n_{i}}$, where
$\mathsf{S}_{i}$ is a diagonal matrix. For a fine element $\tau\subset K_{i}$
which is corresponding to the $s$-th entry of $\mathbb{R}^{n_{i}}$, we can
obtain that $[\mathsf{S}_{i}]_{s,s}=\tilde{\kappa}|_{\tau}\lvert\tau\rvert$,
where $\lvert\tau\rvert$ is the area/volume of $\tau$. Similarly, we can also
rewrite $a_{i}(\cdot,\cdot)$ in Eq. 9 as

$\mathsf{a}_{i}(\mathsf{v},\mathsf{w})=\mathsf{v}\cdot\hat{\mathsf{A}}_{i}\mathsf{w}$

for all $\mathsf{v}$ and $\mathsf{w}$ in $\mathbb{R}^{n}_{i}$, where
$\hat{\mathsf{A}}_{i}\in\mathbb{R}^{n_{i}\times n_{i}}$. Let
$\mathsf{C}_{i}\in\mathbb{R}^{n_{i}\times n}$ be the restriction matrix that
restricts the DoF in the whole domain to $K_{i}$. Thanks to the non-
overlapping property of the partition $\\{K_{i}\\}_{i=1}^{N}$, we can show
that $\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\mathsf{C}_{i}=\mathsf{I}_{n}$,
where $\mathsf{I}_{n}$ is the identity matrix of size $n$. We can have the
following lemma, which is clear from the definitions Eqs. 9 and 7 for a $2$D
uniform structured mesh.

######  Lemma 4.8.

The relation
$\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\hat{\mathsf{A}}_{i}\mathsf{C}_{i}\lesssim\mathsf{A}$
holds.

According to the approximation property of eigenspaces, for any $v_{i}\in
W_{h}(K_{i})$, we can find a $w_{i}\in W_{H}^{\text{c}}(K_{i})$ and $r_{i}\in
W_{h}(K_{i})$ such that $v_{i}=w_{i}+r_{i}$ with

$s_{i}(v_{i}-w_{i},v_{i}-w_{i})=s_{i}(r_{i},r_{i})\leq\frac{1}{\lambda_{L_{i}+1}}a_{i}(v_{i},v_{i}),$

where $\lambda_{L_{i}+1}$ is the $(L_{i}+1)$-th smallest eigenvalue of the
eigenproblem Eq. 11. As illustrated in Section 3, by setting a proper
$\tilde{\kappa}$ and choosing $L_{i}$ large enough, the eigenvalue
$\lambda_{L_{i}+1}$ can be bounded below stably w.r.t. contrast ratios. Let
$\mathsf{E}_{i}\in\mathbb{R}^{n_{i}^{m}\times n_{i}}$ be an extension matrix
that extends the DoF in $K_{i}$ to $K_{i}^{m}$. Our analysis requires the
following assumption.

#### Assumption B

Suppose that $\tilde{\kappa}$ is constructed such that
$\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\lesssim\mathsf{S}_{i}$,
and $L_{i}$ is chosen such that $1/\lambda_{L_{i}+1}\leq\epsilon$ for any $i$,
where $\epsilon$ is a predefined positive number.

To find a $\tilde{\kappa}$ that fulfills Assumption B, a direct way is to lump
all absolute values of the off-diagonal entries into the diagonal entry in
each row. Again, we explain the idea with Fig. 4 for a $2$D uniform structured
mesh. Let a test function $v$ belong to $W_{h}(K_{i})$ that corresponds to a
$\mathsf{v}\in\mathbb{R}^{n_{i}}$. Then we can calculate that

$\displaystyle\mathsf{v}\cdot\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\mathsf{v}=$
$\displaystyle\kappa_{e_{1}}\left(\mathtt{v}_{p,q}-\mathtt{v}_{p,q-1}\right)^{2}\frac{h_{y}}{h_{x}}+\kappa_{e_{2}}\left(\mathtt{v}_{p,q}-\mathtt{v}_{p+1,q}\right)^{2}\frac{h_{x}}{h_{y}}+\kappa_{e_{3}}\mathtt{v}_{p,q}^{2}\frac{h_{x}}{h_{y}}+\kappa_{e_{4}}\mathtt{v}_{p,q}^{2}\frac{h_{y}}{h_{x}}$
$\displaystyle+\left\\{\text{other terms}\right\\},$

where

$\displaystyle\kappa_{e_{1}}=\frac{2}{1/\mathsf{k}_{p,q}+1/\mathsf{k}_{p,q-1}},$
$\displaystyle\kappa_{e_{2}}=\frac{2}{1/\mathsf{k}_{p,q}+1/\mathsf{k}_{p+1,q}},\quad\kappa_{e_{3}}=\frac{2}{1/\mathsf{k}_{p,q}+1/\mathsf{k}_{p,q+1}},$
and
$\displaystyle\kappa_{e_{4}}=\frac{2}{1/\mathsf{k}_{p,q}+1/\mathsf{k}_{p-1,q}}.$

Via the inequalities $\left(\mathtt{v}_{p,q}-\mathtt{v}_{p,q-1}\right)^{2}\leq
2(\mathtt{v}_{p,q}^{2}+\mathtt{v}_{p+1,q}^{2})$ and
$\left(\mathtt{v}_{p,q}-\mathtt{v}_{p,q+1}\right)^{2}\leq
2(\mathtt{v}_{p,q}^{2}+\mathtt{v}_{p,q+1}^{2})$, we can see that

$\displaystyle\mathsf{v}\cdot\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\mathsf{v}\leq$
$\displaystyle\left(2\kappa_{e_{1}}\frac{h_{y}}{h_{x}}+2\kappa_{e_{2}}\frac{h_{x}}{h_{y}}+\kappa_{e_{3}}\frac{h_{y}}{h_{x}}+\kappa_{e_{3}}\frac{h_{x}}{h_{y}}\right)\mathtt{v}_{p,q}^{2}$
$\displaystyle+\left\\{\text{other terms}\right\\}.$

Recalling that $\mathsf{S}_{i}$ is a diagonal matrix with
$[\mathsf{S}_{i}]_{s,s}=\tilde{\kappa}|_{\tau}\lvert\tau\rvert$, we can set

(16)
$\tilde{\kappa}|_{\tau}=\left(2\kappa_{e_{1}}\frac{h_{y}}{h_{x}}+2\kappa_{e_{2}}\frac{h_{x}}{h_{y}}+\kappa_{e_{3}}\frac{h_{y}}{h_{x}}+\kappa_{e_{3}}\frac{h_{x}}{h_{y}}\right)\frac{1}{h_{x}h_{y}},$

where $\tau$ is the fine element in the top-left corner of $K_{i}$. The values
of $\tilde{\kappa}$ on other fine elements could also be found in this way.

$\mathtt{k}_{p,q}$$\mathtt{v}_{p,q}$$\mathtt{k}_{p-1,q}$$\mathtt{k}_{p+1,q}$$\mathtt{v}_{p+1,q}$$\mathtt{k}_{p,q-1}$$\mathtt{v}_{p,q-1}$$\mathtt{k}_{p,q+1}$$K_{i}^{2}$$K_{i}$$e_{1}$$e_{2}$$e_{3}$$e_{4}$
Figure 4: An illustration for a construction of $\tilde{\kappa}$, where
$\mathtt{k}_{p,q}$, $\mathtt{k}_{p-1,q}$, $\mathtt{k}_{p+1,q}$,
$\mathtt{k}_{p,q-1}$ and $\mathtt{k}_{p,q+1}$ are values of $\kappa$ on the
located fine elements respectively, the fine element in the top-left corner of
$K_{i}$ is highlighted with four edges $e_{1}\sim e_{4}$ surrounded,
$\mathtt{v}_{p,q}$, $\mathtt{v}_{p+1,q}$ and $\mathtt{v}_{p,q-1}$ are values
of $v\in W_{h}(K_{i})$ on the located fine elements respectively.

Certainly, if $\tilde{\kappa}$ satisfies Assumption B, then $\tilde{\kappa}$
is positive on each fine element, which implies $\mathsf{S}_{i}$ is
invertible. We can use this information to depict the kernel of
$\mathsf{A}_{0}$ as the following lemma.

######  Lemma 4.9.

Under Assumption B, the dimension of $\ker(\mathsf{A}_{0})$ is $1$, and
$\mathcal{R}$ maps $\ker(\underline{\mathsf{B}})$ into $\ker(\mathsf{A})$,
where $\mathcal{R}$ and $\underline{\mathsf{B}}$ are defined in Eqs. 15 and
14.

######  Proof 4.10.

Let $\mathsf{w}\in\mathbb{R}^{N^{\text{c}}}$ with
$\mathsf{A}_{0}\mathsf{w}=\mathsf{0}$. Then,

$\mathsf{w}\cdot\mathsf{A}_{0}\mathsf{w}=0\Rightarrow\mathsf{w}\cdot\mathsf{R}_{0}\mathsf{A}\mathsf{R}^{\intercal}_{0}\mathsf{w}=0\Rightarrow\sqrt{\mathsf{A}}\mathsf{R}^{\intercal}_{0}\mathsf{w}=\mathsf{0}\Rightarrow\mathsf{R}^{\intercal}_{0}\mathsf{w}\in\ker(\mathsf{A}).$

This gives us that the function $w_{h}\in W_{h}$ represented by
$\mathsf{R}^{\intercal}_{0}\mathsf{w}$ is constant in $\Omega$. On each coarse
element $K_{i}$, the smallest eigenvalue is $0$ with a corresponding
eigenvector is a constant vector. Due to $\mathsf{S}_{i}$ is invertible, the
eigenspace corresponding to $0$ has a dimension of one. Note that
$\\{K_{i}\\}_{i=1}^{N}$ is non-overlapping, then restriction of $w_{h}$ on
$K_{i}$ has only one DoF, which hence determines $\mathsf{w}$. Although that
$\mathsf{w}$ may not be a constant vector, we have
$\dim(\ker(\mathsf{A}_{0}))=1$. It is clear that $\mathcal{R}$ Eq. 15 maps
$\ker(\underline{\mathsf{B}})$ into $\ker(\mathsf{A})$ because
$\mathsf{\mathsf{A}}_{i}$ is invertible and $\ker(\mathsf{A}_{0})$ has been
fully understood.

The following lemma provides an estimation of $C_{T}$ in Eq. 12.

######  Lemma 4.11.

Under Assumption A and B, for all
$\mathsf{u}\in\operatorname{im}(\mathsf{A})$, there exists a
$\underline{\mathsf{v}}\in\operatorname{im}(\underline{\mathsf{B}})$ with
$\mathsf{u}=\mathcal{R}\underline{\mathsf{v}}$ and

$\frac{1}{2+(8C_{\text{os}}+1)\epsilon}\,\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}})\leq\mathsf{a}(\mathsf{u},\mathsf{u}),$

where $\underline{\mathsf{B}}$, $\mathcal{R}$ and $\mathsf{b}(\cdot,\cdot)$
are defined in Eqs. 15 and 14.

######  Proof 4.12.

We decompose $\mathsf{C}_{i}\mathsf{u}$ as
$\mathsf{C}_{i}\mathsf{u}=\mathsf{w}_{i}+\mathsf{r}_{i}$, where the function
representation of $\mathsf{w}_{i}$ belongs to $W_{H}^{\text{c}}(K_{i})$.
According to the approximation property of eigenspaces, we can find a
$\mathsf{w}_{i}\in\mathbb{R}^{n_{i}}$ such that
$\mathsf{s}_{i}(\mathsf{r}_{i},\mathsf{r}_{i})$ is minimized with an
estimation

(17)
$\mathsf{s}_{i}(\mathsf{r}_{i},\mathsf{r}_{i})=\mathsf{r}_{i}\cdot\mathsf{S}_{i}\mathsf{r}_{i}\leq\epsilon\mathsf{a}_{i}(\mathsf{C}_{i}\mathsf{u},\mathsf{C}_{i}\mathsf{u})=\epsilon\mathsf{u}\cdot\mathsf{C}_{i}^{\intercal}\hat{\mathsf{A}}_{i}\mathsf{C}_{i}\mathsf{u}.$

Recalling that
$\mathsf{u}=\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\mathsf{C}_{i}\mathsf{u}$,
we have

$\mathsf{u}=\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\left(\mathsf{w}_{i}+\mathsf{r}_{i}\right)=\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\mathsf{w}_{i}+\sum_{i=1}^{N}\mathsf{R}_{i}^{\intercal}\mathsf{r}_{i}^{\prime},$

where $\mathsf{r}_{i}^{\prime}=\mathsf{E}_{i}\mathsf{r}_{i}$. Moreover, from
constructions of the global coarse space $W_{H}^{\text{c}}$ and
$\mathsf{R}_{0}$, there exists a $\mathsf{w}\in\mathbb{R}^{N^{\text{c}}}$ such
that
$\sum_{i=1}^{N}\mathsf{C}_{i}^{\intercal}\mathsf{w}_{i}=\mathsf{R}_{0}^{\intercal}\mathsf{w}$.
We hence derive a vector $\underline{\mathsf{v}}\in\underline{V}$ as

$\underline{\mathsf{v}}=[\mathsf{w},\mathsf{r}_{1}^{\prime},\dots,\mathsf{r}_{N}^{\prime}]$

that satisfies $\mathsf{u}=\mathcal{R}\underline{\mathsf{v}}$. Thanks to Lemma
4.9, it is easy to check that if $\mathsf{u}\in\operatorname{im}(\mathsf{A})$,
the constructed vector $\underline{\mathsf{v}}$ belongs to
$\operatorname{im}(\underline{\mathsf{B}})$.

Note that

$\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}})=\mathsf{w}\cdot\mathsf{A}_{0}\mathsf{w}+\sum_{i=1}^{N}\mathsf{r}_{i}^{\prime}\cdot\mathsf{A}_{i}\mathsf{r}_{i}^{\prime}=\mathsf{a}(\mathsf{R}_{0}^{\intercal}\mathsf{w},\mathsf{R}_{0}^{\intercal}\mathsf{w})+\sum_{i=1}^{N}\mathsf{r}_{i}\cdot\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\mathsf{r}_{i}.$

We first need to bound
$\mathsf{a}(\mathsf{R}_{0}^{\intercal}\mathsf{w},\mathsf{R}_{0}^{\intercal}\mathsf{w})$.
Employing techniques in the proof of Lemma 4.6, we have

$\displaystyle\mathsf{a}(\mathsf{R}_{0}^{\intercal}\mathsf{w},\mathsf{R}_{0}^{\intercal}\mathsf{w})\leq$
$\displaystyle
2\left\\{\mathsf{a}(\mathsf{u},\mathsf{u})+\mathsf{a}(\sum_{i}^{N}\mathsf{R}_{i}^{\intercal}\mathsf{r}_{i}^{\prime},\sum_{i}^{N}\mathsf{R}_{i}^{\intercal}\mathsf{r}_{i}^{\prime})\right\\}$
$\displaystyle\leq$ $\displaystyle
2\left\\{\mathsf{a}(\mathsf{u},\mathsf{u})+4C_{\text{os}}\sum_{i=1}^{N}\mathsf{r}_{i}^{\prime}\cdot\mathsf{A}_{i}\mathsf{r}_{i}^{\prime}\right\\}$
$\displaystyle=$ $\displaystyle
2\left\\{\mathsf{a}(\mathsf{u},\mathsf{u})+4C_{\text{os}}\sum_{i=1}^{N}\mathsf{r}_{i}\cdot\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\mathsf{r}_{i}\right\\}.$

Combining Assumption B with Eq. 17, we arrive at

$\displaystyle\sum_{i=1}^{N}\mathsf{r}_{i}\cdot\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\mathsf{r}_{i}\leq$
$\displaystyle\sum_{i=1}^{N}\mathsf{r}_{i}\cdot\mathsf{S}_{i}\mathsf{r}_{i}\leq\epsilon\sum_{i=1}^{N}\mathsf{u}\cdot\mathsf{C}_{i}^{\intercal}\hat{\mathsf{A}}_{i}\mathsf{C}_{i}\mathsf{u}$
$\displaystyle\leq$
$\displaystyle\epsilon\mathsf{u}\cdot\mathsf{A}\mathsf{u}=\epsilon\mathsf{a}(\mathsf{u},\mathsf{u}),$

where the last line is from Lemma 4.8. Therefore, we derive that

$\mathsf{b}(\underline{\mathsf{v}},\underline{\mathsf{v}})\leq
2\left\\{1+4C_{\text{os}}\epsilon\right\\}\mathsf{a}(\mathsf{u},\mathsf{u})+\epsilon\mathsf{a}(\mathsf{u},\mathsf{u})\leq\left\\{2+(8C_{\text{os}}+1)\epsilon\right\\}\mathsf{a}(\mathsf{u},\mathsf{u}),$

which leads the completion of the proof.

Combining Lemmas 4.6, 4.11, and 4.9, we reach the main theorem of this
section.

######  Theorem 4.13.

Under Assumption A and B, the inequality

$\operatorname{cond}(\mathsf{P}^{-1}\mathsf{A})\leq\left(2+(8C_{\text{os}}+1)\epsilon\right)\left(1+4C_{\text{os}}\right)$

holds.

##  5 Numerical experiments

We implement our method with PETSc [8]. We set the unit cube $(0,1)^{3}$ as
the computational domain with uniform structured meshes for all experiments
(expect for the SPE10 problem). The domain decomposition is handled by PETSc’s
DMDA module, i.e., each MPI process owns a non-overlapping rectangular cuboid
region plus $m$ layers of neighboring ghost points for inter-process
communications. Meanwhile, each cuboid is further divided into $\mathtt{sd}$
parts in $x$-, $y$\- and $z$-direction, which gives a non-overlapping
decomposition of the whole domain into total
$\mathtt{sd}^{3}\times\mathtt{proc}$ subdomains (also called coarse elements),
where $\mathtt{proc}$ stands for the number of MPI processes, and this also
signifies $N=\mathtt{sd}^{3}\times\mathtt{proc}$. Note that in our computing
settings, coarse elements are not of the same size, whilst $H$ the minimal
size could be roughly estimated as
$1.0/(\mathtt{sd}\times\sqrt[3]{\mathtt{proc}})$. In each coarse element, we
solve the eigenvalue problem Eq. 11 with SLEPc [32] and obtain $L_{\star}$
eigenvectors. For local and global solvers, we take CHOLMOD [14] to factorize
$\mathsf{A}_{i}$ and MUMPS [3] to factorize $\mathsf{A}_{0}$.

We place $4$ long singular sources at $4$ corner points on $x\times y$-plane
and an opposite sign long singular source at the middle of $x\times y$-plane
as the source term $f$ of the model Eq. 2, which mimics the well condition in
the SPE10 model (ref. [16]). We set the default parameters for the GMRES
solver provided by PETSc to solve the linear algebraic system Eq. 5, and the
convergence criterion is that $l^{2}$-norm of the residual reaches $10^{-5}$.
All computational elapsed time records in numerical experiments are measured
from an HPC cluster that each computing node is configured with dual Intel®
Xeon® Gold 6140 CPUs (36 cores) and 192GB memory. In the following reports, We
denote by $\mathtt{iter}$ the number of GMRES iterations, $\mathtt{DoF}$ the
number of total DoF.

We address several different settings in numerical experiments compared with
Assumption B in Section 4: first, we preset a fixed number $L_{\star}$ of
eigenvectors for all eigenproblems; second, we simply choose
$\tilde{\kappa}=\kappa$ rather than the sophisticated construction in Eq. 16.
We explain those modifications here. An important topic in designing parallel
algorithms is the so-called load balancing. In our situation, Theorem 4.13
implies that determining $L_{i}$ the number of eigenvectors adaptively
according to the threshold for eigenvalues may achieve a theoretical bound for
the condition number. However, for some challenging permeability profiles, the
variance of $\\{L_{i}\\}_{i=1}^{N}$ may be considerably large, which causes a
load imbalance among processes and will tremendously degrade the overall
performance. Moreover, thanks to the non-overlapping property, we could
directly construct $L_{\star}$ global vectors defined on the fine mesh for
storing all eigenvectors (coarse bases). Then, we can directly leverage the
functionality from DMDA for handling communications in building the global
solver $\mathsf{A}_{0}$, which greatly facilitates our implementation. For
choosing $\tilde{\kappa}=\kappa$, recalling in Eq. 16, we obtain an expression
for $\tilde{\kappa}|_{\tau}$ with $\kappa_{e_{1}}$, $\kappa_{e_{2}}$,
$\kappa_{e_{3}}$ and $\kappa_{e_{4}}$, and we can show that $\kappa_{e_{1}}$,
$\kappa_{e_{2}}$, $\kappa_{e_{3}}$ and $\kappa_{e_{4}}$ can all be bounded by
$2\mathtt{k}_{p,q}$ which is $2\kappa|_{\tau}$. This observation implies that
by taking $\tilde{\kappa}=C\kappa$ with a modest large constant $C$, the
inequality
$\mathsf{E}_{i}^{\intercal}\mathsf{A}_{i}\mathsf{E}_{i}\lesssim\mathsf{S}_{i}$
is obtainable. Since we only require a fixed number of eigenvectors, the
constant $C$ does not affect the resulting eigenspace (also called the local
coarse space). Again, taking $\tilde{\kappa}=\kappa$ is easier to implement in
contrast to Eq. 16, because the latter needs to cook up different formulas
depending on whether the fine element touches $\partial K_{i}$.

###  5.1 Contrast robustness

In this subsection, we will test our preconditioner with two types of media,
that is, one contains long channels and the other contains fractures. We will
fix the coefficient $\kappa=1$ in the background region, while increasing
values of $\kappa$ in channels and fractures to examine robustness w.r.t.
contrast ratios. All models have the same resolution as $\mathtt{DoF}=512^{3}$
and are solved with the same number of MPI processes as $\mathtt{proc}=8^{3}$.

In the first type of media, we place long channels that penetrate the whole
domain in $y$-direction, that is $\kappa(x,y,z)=\kappa(x,0,z)$ for all $y$.
The resolution of $x\times z$-plane is of $512\times 512$, and we arrange
periodic cells that are in a size of $16\times 16$ along $x$\- and
$z$-direction. The permeability field $\kappa$ is determined by two parameters
$\alpha$ and $\mathtt{cr}$ with $\kappa=10^{\alpha\times\mathtt{cr}}$, where
distributions of $\alpha$ in a periodic cell are demonstrated in the left part
of Fig. 5. We consider four different configurations—“2-ch. cfg.”, “3-ch.
cfg.”, “4-ch. cfg.” and “5-ch. cfg.”—which are named by the number of channels
inside a periodic cell. The number of GMRES iterations $\mathtt{iter}$ against
$L_{\star}$ under different $\mathtt{cr}$ is plotted in the right part of Fig.
5, and note that we set $m=3$ and $\mathtt{sd}=4$ for those tests. From the
graphs, we could see that $\mathtt{iter}$ can remain robust by enriching
coarse spaces. For example, for each configuration, when $L_{\star}=6$ or $7$,
all values of $\mathtt{iter}$ are all close with contrast ratios ranging from
$10^{3}\sim 10^{6}$. Moreover, we could observe that problems are tougher with
more channels (e.g., for “4-ch. cfg.” and “5-ch. cfg.” with $\mathtt{cr}=12$,
the solver cannot finish within $1000$ iterations). Again, our preconditioner
can achieve robust performance by incorporating a modest number of
eigenvectors ($L_{\star}=6$) into coarse spaces.

Figure 5: (left part) the distributions of $\alpha$ for different
configurations in a periodic cell, (right part) The number of GMRES iterations
$\mathtt{iter}$ against $L_{\star}$ under different $\mathtt{cr}$.

The second type of media are visualized in the left part of Fig. 6 that
contains two configurations, and the original models are in a resolution of
$256^{3}$ while we periodically duplicate them into $512^{3}$ the resolution
of computational domains. Note that in Fig. 6, blue regions refer to
fractures, where we will place significantly higher permeability as
$\kappa=10^{\mathtt{cr}}$ in to reflect that flows move fast in fractures, and
we simply fixed $\kappa=1$ in red regions. We emphasized that those models
are, indeed, challenging to solve. To support this assertion, we adopt
“gamg”—the default algebraic multigrid preconditioner by PETSc—to GMRES and
record $\mathtt{iter}$ and elapsed times in Table 1. It is found that the
performance of “gamg” gradually deteriorates as increasing $\mathtt{cr}$, and
there are errors of breakdown appearing when $\mathtt{cr}=6$ and $8$. However,
our preconditioner exhibits robustness w.r.t. $\mathtt{cr}$ from the data in
Table 1 with a set of parameters $(L_{\star},m,\mathtt{sd})=(4,2,4)$. The
right part of Fig. 6 shows $\mathtt{iter}$ against $L_{\star}$ under different
$\mathtt{cr}$ for two fracture configurations with $m=3$ and $\mathtt{sd}=4$,
where the contrast ratios in those tests range from $1$ to $10^{10}$. The
plots again illustrate that the number of iterations decreases as expected
with increasing the dimension of coarse spaces, and a modest large $L_{\star}$
is necessary to achieve robustness for those high-contrast problems.

Figure 6: (left part) two binary-valued data sets of a size of $256^{3}$, where the regions of blue color represent fractures, and the whole computational domain $512^{3}$ is constituted by periodically duplicating a $256^{3}$ configuration; (right part) plots of $\mathtt{iter}$ w.r.t. different $\mathtt{en}$ under different $\mathtt{cr}$. Table 1: Iteration numbers and elapsed times (seconds) in a format $\mathtt{iter}(\text{time})$ w.r.t. different $\mathtt{cr}$, where “gamg” stands for PETSc’s default algebraic multigrid preconditioner and “$-$” for breakdown error occurring, and the coefficient profile used is from “B-frac. cfg.” in Fig. 6. $\mathtt{cr}$ | $0$ | $2$ | $4$ | $6$ | $8$
---|---|---|---|---|---
gamg | $9(27.8)$ | $18(33.1)$ | $99(77.8)$ | $-$ | $-$
$(L_{\star},m,\mathtt{sd})=(4,2,4)$ | $23(40.8)$ | $44(61.1)$ | $55(77.2)$ | $60(81.5)$ | $55(86.8)$

###  5.2 Scalability

In this subsection, we fix $\mathtt{sd}=4$, $m=3$ and $L_{\star}=4$ while
change $\mathtt{DoF}$ and $\mathtt{proc}$. Moreover, the permeability fields
are generated from the periodical cell “2-ch. cfg.” in Fig. 5 with
$\mathtt{cr}=6$. Compared to Algorithm 1, we divide the elapsed time of
obtaining solutions into three parts: Pre-0 (construct and factorize
$\mathsf{A}_{i}$, solve local eigenvalue problems), Pre-1 (construct and
factorize $\mathsf{A}_{0}$) and Ite (the iteration phase). The results of
strong and weak scalability tests are presented in the left and right part of
Fig. 7 respectively. For examining strong scalability, we set
$\mathtt{DoF}=512^{3}$ and increase $\mathtt{proc}$ from $6^{3}$ to $10^{3}$,
and a three-fold acceleration can be observed. Recall that our domain
decomposition is related to the number of processes and $H\sim
1.0/(\mathtt{sd}\times\sqrt[3]{\mathtt{proc}})$, we can also notice a decrease
in $\mathtt{iter}$ with a large $\mathtt{proc}$, and this phenomenon agrees
with traditional theories that a smaller $H$ can result in a better condition
number. To test weak scalability, we fix the ratio between $\mathtt{DoF}$ and
$\mathtt{proc}$, and the performance of $\mathtt{DoF}=640^{3}$ degrades to
$70\%$ of $\mathtt{DoF}=256^{3}$, while $\mathtt{iter}$ remains unchanged. The
bottleneck of our method for obtaining better scalability performance is
global solvers, or more precisely, factorize $\mathsf{A}_{0}$. We could see
this argument more clearly from the elapsed times in the Pre-1 rows in Fig. 7,
where both cases show a dramatic increase with more MPI processes. Tuning
direct solvers for sparse linear systems is a well-known difficulty [48], and
a remedy (also in our future work) could be replacing the factorization of
$\mathsf{A}_{0}$ with an in-exact but scalable solver.

Figure 7: (left part) elapsed times and iteration numbers in strong
scalability tests with $\mathtt{DoF}=512^{3}$, (right part) elapsed times and
iteration numbers in weak scalability tests with a fixed ratio between
$\mathtt{DoF}$ and $\mathtt{proc}$.

###  5.3 Parameter tests

In this subsection, we fix $\mathtt{DoF}=512^{3}$, $\mathtt{proc}=8^{3}$ to
study the influence of $\mathtt{sd}$, $m$ and $L_{\star}$ on the solver. The
permeability field $\kappa$ is constructed from “B-frac. cfg.” in Fig. 6 with
$\mathtt{cr}=6$. Although Theorem 4.13 does not depict how the number of
oversampling layers $m$ affects the condition number, Table 2 indeed confirms
that $\mathtt{iter}$ could be smaller with a larger $m$. However, as we have
mentioned before, a larger $m$ means more communications and will degrade the
final performance, which could also be validated from the elapsed times in
Table 2 that the best time record occurs with $m=1$ and $L_{\star}=4$.

Table 2: Iteration numbers and elapsed times (seconds) in a format $\mathtt{iter}(\text{time})$ w.r.t. different $m$ and $L_{\star}$ with $\mathtt{sd}=4$.  $L_{\star}$ | $3$ | $4$ | $5$ | $6$ | $7$
---|---|---|---|---|---
$m=1$ | $107(89.6)$ | $83(80.9)$ | $80(87.3)$ | $73(93.7)$ | $58(107.5)$
$m=2$ | $78(94.7)$ | $60(81.5)$ | $54(100.0)$ | $50(86.7)$ | $47(91.3)$
$m=3$ | $68(120.9)$ | $53(105.6)$ | $47(102.6)$ | $44(103.5)$ | $42(110.3)$
$m=4$ | $56(160.6)$ | $48(138.2)$ | $42(130.8)$ | $39(134.8)$ | $37(134.1)$

Note that in our settings, we introduce $\mathtt{sd}$ to further divide the
computational domain owned by each process. If we just set $\mathtt{sd}=1$ for
restricting the time consumed on local solvers, we are left to use more MPI
processes, which usually is unaffordable due to limited computational
resources. The results of experiments on $\mathtt{sd}$ are presented in Table
3. Just as the strong scalability test in Fig. 7, a larger $\mathtt{sd}$ means
a smaller $H$, which decreases the number of iterations. However, a larger
$\mathtt{sd}$ can increase the size of $\mathsf{A}_{0}$, and the final
performance will be burdened with global solvers. The key observation from
Table 3 is that to efficiently implement two-level preconditioners, balancing
the efforts on local and global solvers is crucial.

Table 3: Iteration numbers and elapsed times (seconds) in a format $\mathtt{iter}(\text{time})$ w.r.t. different $\mathtt{sd}$ and $L_{\star}$ with $m=3$. $L_{\star}$ | $3$ | $4$ | $5$ | $6$ | $7$
---|---|---|---|---|---
$\mathtt{sd}=1$ | $-$ | $173(495.4)$ | $153(510.9)$ | $138(442.0)$ | $133(429.0)$
$\mathtt{sd}=2$ | $-$ | $130(186.6)$ | $127(189.4)$ | $68(123.1)$ | $64(121.1)$
$\mathtt{sd}=4$ | $68(120.9)$ | $53(105.6)$ | $47(102.6)$ | $44(103.5)$ | $42(110.3)$
$\mathtt{sd}=8$ | $28(185.0)$ | $24(202.1)$ | $22(268.1)$ | $22(273.1)$ | $22(331.3)$

Next, we consider the SPE10 problem [16], which is a widely used benchmark
model for testing various model reduction methods, and the adopted
permeability field is visualized in the left part of Fig. 8. This model
exhibits strong heterogeneity with a magnitude of $10^{7}$. The computed
pressure field is demonstrated in the right part of Fig. 8. The $\mathtt{DoF}$
in this model is around $1$ million, which is relatively small in comparison
with previous experiments, and we just report several results with different
$L_{\star}$ in Table 4. We could observe again that a larger $L_{\star}$ can
bring a smaller $\mathtt{iter}$, while the elapsed times may get improved by
exploring other settings of $\mathtt{proc}$, $m$ and $\mathtt{sd}$.

Figure 8: (left part) the coefficient profile of the SPE10 model, (right part) The pressure distribution. Table 4: Iteration numbers and elapsed times (seconds) in a format $\mathtt{iter}(\text{time})$ w.r.t. different $L_{\star}$ for the SPE10 model with $\mathtt{proc}=56$, $m=3$ and $\mathtt{sd}=2$. $L_{\star}$ |  $1$ |  $2$ |  $3$ |  $4$ |  $5$ |  $6$ |  $7$
---|---|---|---|---|---|---|---
$\mathtt{iter}(\text{wall time})$ |  $-$ |  $-$ |  $-$ |  $49(6.5)$ |  $39(5.4)$ |  $33(4.8)$ |  $29(4.3)$

##  6 Conclusion

In this paper, we developed a two-level overlapping preconditioner for
accelerating iterative solvers for linear algebraic systems derived from
discretizing Darcy flow. To prepare for practical two-phase flow simulations,
the discretization that we focus on is based on the mixed form that is
superior in mass conservation compared with the common second-order form, and
the no-flux boundary condition is also considered rather than for example, the
Dirichlet boundary condition. Inheriting the methodology of GMsFEMs, we solve
general eigenvalue problems on subdomains to obtain leading eigenvectors for
constructing the coarse space. Our construction distinguishes itself from
existing literature by working on non-overlapping coarse elements, which to
the largest extent saves communications and guarantees parallel efficiency.
Our manner of building local solvers is not a direct copy of additive Schwarz
methods on Galerkin discretization, and actually, it has been explained to be
suitable for parallel implementations. To handle the singularity from the no-
flux boundary condition, we tweaked the fictitious space lemma and obtained an
a priori estimation for the conditioner number of the preconditioner under
several assumptions. Our implementation is facilitated by PETSc, and numerical
reports support the claim that the proposed preconditioner is robust with high
contrast coefficient profiles. Furthermore, the strong and weak scalability
tests inspire a possible improvement—replacing the factorization in the global
solver with an in-exact but scalable one–that will be provided in our future
work.

## References

  * [1] J. r. E. Aarnes, On the use of a mixed multiscale finite element method for greater flexibility and increased speed or improved accuracy in reservoir simulation, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 2 (2004), pp. 421–439, https://doi.org/10.1137/030600655.
  * [2] P. R. Amestoy, A. Buttari, J.-Y. L’Excellent, and T. Mary, Performance and scalability of the block low-rank multifrontal factorization on multicore architectures, ACM Transactions on Mathematical Software, 45 (2019), https://doi.org/10.1145/3242094.
  * [3] P. R. Amestoy, I. S. Duff, J.-Y. L'Excellent, and J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 15–41, https://doi.org/10.1137/s0895479899358194.
  * [4] T. Arbogast, G. Pencheva, M. F. Wheeler, and I. Yotov, A multiscale mortar mixed finite element method, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 6 (2007), pp. 319–346, https://doi.org/10.1137/060662587.
  * [5] T. Arbogast, M. F. Wheeler, and I. Yotov, Mixed finite elements for elliptic problems with tensor coefficients as cell-centered finite differences, SIAM Journal on Numerical Analysis, 34 (1997), pp. 828–852, https://doi.org/10.1137/S0036142994262585.
  * [6] T. Arbogast and H. Xiao, A multiscale mortar mixed space based on homogenization for heterogeneous elliptic problems, SIAM Journal on Numerical Analysis, 51 (2013), pp. 377–399, https://doi.org/10.1137/120874928.
  * [7] F. Auricchio, L. B. a. da Veiga, F. Brezzi, and C. Lovadina, Mixed finite element methods, John Wiley & Sons, Ltd, 2017, pp. 1–53, https://doi.org/https://doi.org/10.1002/9781119176817.ecm2004.
  * [8] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, PETSc/TAO users manual, Tech. Report ANL-21/39 - Revision 3.17, Argonne National Laboratory, 2022\.
  * [9] D. Boffi, F. Brezzi, and M. Fortin, Mixed finite element methods and applications, vol. 44 of Springer Series in Computational Mathematics, Springer, Heidelberg, 2013, https://doi.org/10.1007/978-3-642-36519-5.
  * [10] W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second ed., 2000, https://doi.org/10.1137/1.9780898719505.
  * [11] P. N. Brown and H. F. Walker, GMRES on (nearly) singular systems, SIAM Journal on Matrix Analysis and Applications, 18 (1997), pp. 37–51, https://doi.org/10.1137/S0895479894262339.
  * [12] J. G. Calvo and O. B. Widlund, An adaptive choice of primal constraints for BDDC domain decomposition algorithms, Electronic Transactions on Numerical Analysis, 45 (2016), pp. 524–544.
  * [13] J. Chen, E. T. Chung, Z. He, and S. Sun, Generalized multiscale approximation of mixed finite elements with velocity elimination for subsurface flow, Journal of Computational Physics, 404 (2020), pp. 109133, 23, https://doi.org/10.1016/j.jcp.2019.109133.
  * [14] Y. Chen, T. A. Davis, W. W. Hager, and S. Rajamanickam, Algorithm 887: CHOLMOD, supernodal sparse Cholesky factorization and update/downdate, Association for Computing Machinery. Transactions on Mathematical Software, 35 (2008), https://doi.org/10.1145/1391989.1391995.
  * [15] Z. Chen and T. Y. Hou, A mixed multiscale finite element method for elliptic problems with oscillating coefficients, Mathematics of Computation, 72 (2003), pp. 541–576, https://doi.org/10.1090/S0025-5718-02-01441-2.
  * [16] M. A. Christie and M. J. Blunt, Tenth SPE comparative solution project: a comparison of upscaling techniques, SPE Reservoir Evaluation & Engineering, 4 (2001), pp. 308–317, https://doi.org/10.2118/72469-PA, https://doi.org/10.2118/72469-PA.
  * [17] E. T. Chung, Y. Efendiev, and C. S. Lee, Mixed generalized multiscale finite element methods and applications, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 13 (2015), pp. 338–366, https://doi.org/10.1137/140970574.
  * [18] V. Dolean, P. Jolivet, and F. Nataf, An introduction to domain decomposition methods, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015, https://doi.org/10.1137/1.9781611974065.ch1.  Algorithms, theory, and parallel implementation.
  * [19] V. Dolean, F. Nataf, R. Scheichl, and N. Spillane, Analysis of a two-level Schwarz method with coarse spaces based on local Dirichlet-to-Neumann maps, Computational Methods in Applied Mathematics, 12 (2012), pp. 391–414, https://doi.org/10.2478/cmam-2012-0027.
  * [20] L. J. Durlofsky, Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resources Research, 27 (1991), pp. 699–708, https://doi.org/https://doi.org/10.1029/91WR00107.
  * [21] Y. Efendiev, J. Galvis, and T. Y. Hou, Generalized multiscale finite element methods (GMsFEM), Journal of Computational Physics, 251 (2013), pp. 116–135, https://doi.org/10.1016/j.jcp.2013.04.045.
  * [22] Y. Efendiev, J. Galvis, R. Lazarov, and J. Willems, Robust domain decomposition preconditioners for abstract symmetric positive definite bilinear forms, ESAIM. Mathematical Modelling and Numerical Analysis, 46 (2012), pp. 1175–1199, https://doi.org/10.1051/m2an/2011073.
  * [23] J. Galvis and Y. Efendiev, Domain decomposition preconditioners for multiscale flows in high-contrast media, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 8 (2010), pp. 1461–1483, https://doi.org/10.1137/090751190.
  * [24] J. Galvis and Y. Efendiev, Domain decomposition preconditioners for multiscale flows in high contrast media: reduced dimension coarse spaces, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 8 (2010), pp. 1621–1644, https://doi.org/10.1137/100790112.
  * [25] I. G. Graham and R. Scheichl, Robust domain decomposition algorithms for multiscale PDEs, Numerical Methods for Partial Differential Equations. An International Journal, 23 (2007), pp. 859–878, https://doi.org/10.1002/num.20254.
  * [26] M. Griebel and P. Oswald, On the abstract theory of additive and multiplicative Schwarz algorithms, Numerische Mathematik, 70 (1995), pp. 163–180, https://doi.org/10.1007/s002110050115.
  * [27] H. Hajibeygi and P. Jenny, Multiscale finite-volume method for parabolic problems arising from compressible multiphase flow in porous media, Journal of Computational Physics, 228 (2009), pp. 5129–5147, https://doi.org/10.1016/j.jcp.2009.04.017.
  * [28] Z. He, H. Chen, J. Chen, and Z. Chen, Generalized multiscale approximation of a mixed finite element method with velocity elimination for Darcy flow in fractured porous media, Computer Methods in Applied Mechanics and Engineering, 381 (2021), p. 113846, https://doi.org/10.1016/j.cma.2021.113846.
  * [29] Z. He, E. T. Chung, J. Chen, and Z. Chen, Adaptive generalized multiscale approximation of a mixed finite element method with velocity elimination, Computational Geosciences. Modeling, Simulation and Data Analysis, 25 (2021), pp. 1681–1708, https://doi.org/10.1007/s10596-021-10068-9.
  * [30] A. Heinlein, A. Klawonn, J. Knepper, and O. Rheinbach, Multiscale coarse spaces for overlapping Schwarz methods based on the ACMS space in 2D, Electronic Transactions on Numerical Analysis, 48 (2018), pp. 156–182, https://doi.org/10.1553/etna_{v}{o}{l}48s156.
  * [31] A. Heinlein, A. Klawonn, J. Knepper, and O. Rheinbach, Adaptive GDSW coarse spaces for overlapping Schwarz methods in three dimensions, SIAM Journal on Scientific Computing, 41 (2019), pp. A3045–A3072, https://doi.org/10.1137/18M1220613.
  * [32] V. Hernandez, J. E. Roman, and V. Vidal, SLEPc, ACM Transactions on Mathematical Software, 31 (2005), pp. 351–362, https://doi.org/10.1145/1089014.1089019.
  * [33] T. Y. Hou and X.-H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, Journal of Computational Physics, 134 (1997), pp. 169–189, https://doi.org/10.1006/jcph.1997.5682.
  * [34] H. H. Kim, E. Chung, and J. Wang, BDDC and FETI-DP preconditioners with adaptive coarse spaces for three-dimensional elliptic problems with oscillatory and high contrast coefficients, Journal of Computational Physics, 349 (2017), pp. 191–214, https://doi.org/10.1016/j.jcp.2017.08.003.
  * [35] H. H. Kim, E. Chung, and J. Wang, BDDC and FETI-DP algorithms with a change of basis formulation on adaptive primal constraints, Electronic Transactions on Numerical Analysis, 49 (2018), pp. 64–80, https://doi.org/10.1553/etna_vol49s64.
  * [36] A. Klawonn, P. Radtke, and O. Rheinbach, FETI-DP methods with an adaptive coarse space, SIAM Journal on Numerical Analysis, 53 (2015), pp. 297–320, https://doi.org/10.1137/130939675.
  * [37] A. Klawonn, P. Radtke, and O. Rheinbach, A comparison of adaptive coarse spaces for iterative substructuring in two dimensions, Electronic Transactions on Numerical Analysis, 45 (2016), pp. 75–106.
  * [38] R. Li, H. Yang, and C. Yang, Parallel multilevel restricted Schwarz preconditioners for implicit simulation of subsurface flows with Peng-Robinson equation of state, Journal of Computational Physics, 422 (2020), p. 109745, https://doi.org/10.1016/j.jcp.2020.109745.
  * [39] I. Lunati and P. Jenny, Multiscale finite-volume method for compressible multiphase flow in porous media, Journal of Computational Physics, 216 (2006), pp. 616–636, https://doi.org/10.1016/j.jcp.2006.01.001.
  * [40] I. Lunati and P. Jenny, Treating highly anisotropic subsurface flow with the multiscale finite-volume method, Multiscale Modeling & Simulation. A SIAM Interdisciplinary Journal, 6 (2007), pp. 308–318, https://doi.org/10.1137/050638928.
  * [41] J. Mandel and B. Sousedík, Adaptive selection of face coarse degrees of freedom in the BDDC and the FETI-DP iterative substructuring methods, Computer Methods in Applied Mechanics and Engineering, 196 (2007), pp. 1389–1399, https://doi.org/10.1016/j.cma.2006.03.010.
  * [42] A. M. Manea, J. Sewall, and H. A. Tchelepi, Parallel multiscale linear solver for highly detailed reservoir models, SPE Journal, 21 (2016), pp. 2062–2078, https://doi.org/10.2118/173259-PA, https://doi.org/10.2118/173259-PA.
  * [43] F. Nataf, H. Xiang, and V. Dolean, A two level domain decomposition preconditioner based on local Dirichlet-to-Neumann maps, Comptes Rendus Mathematique, 348 (2010), pp. 1163–1167, https://doi.org/10.1016/j.crma.2010.10.007.
  * [44] F. Nataf, H. Xiang, V. Dolean, and N. Spillane, A coarse space construction based on local Dirichlet-to-Neumann maps, SIAM Journal on Scientific Computing, 33 (2011), pp. 1623–1642, https://doi.org/10.1137/100796376.
  * [45] Y. Saad, Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, Philadelphia, PA, second ed., 2003, https://doi.org/10.1137/1.9780898718003.
  * [46] M. Sarkis, Nonstandard coarse spaces and Schwarz methods for elliptic problems with discontinuous coefficients using non-conforming elements, Numerische Mathematik, 77 (1997), pp. 383–406, https://doi.org/10.1007/s002110050292.
  * [47] A. Toselli and O. Widlund, Domain decomposition methods—algorithms and theory, vol. 34 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 2005, https://doi.org/10.1007/b137868.
  * [48] S. Wang, M. V. de Hoop, and J. Xia, On 3D modeling of seismic wave propagation via a structured parallel multifrontal direct Helmholtz solver, Geophysical Prospecting, 59 (2011), pp. 857–873, https://doi.org/10.1111/j.1365-2478.2011.00982.x.
  * [49] Y. Wang, H. Hajibeygi, and H. A. Tchelepi, Algebraic multiscale solver for flow in heterogeneous porous media, Journal of Computational Physics, 259 (2014), pp. 284–303, https://doi.org/10.1016/j.jcp.2013.11.024.
  * [50] X. H. Wu, Y. Efendiev, and T. Y. Hou, Analysis of upscaling absolute permeability, Discrete and Continuous Dynamical Systems. Series B. A Journal Bridging Mathematics and Sciences, 2 (2002), pp. 185–204, https://doi.org/10.3934/dcdsb.2002.2.185.
  * [51] H. Yang, S. Sun, Y. Li, and C. Yang, A fully implicit constraint-preserving simulator for the black oil model of petroleum reservoirs, Journal of Computational Physics, 396 (2019), pp. 347–363, https://doi.org/10.1016/j.jcp.2019.05.038.
  * [52] H. Yang, Z. Zhu, and J. Kou, A minimum-type nonlinear complementarity simulator with constrained pressure residual (CPR) methods for wormhole propagation in carbonate acidization, Journal of Computational Physics, 473 (2023), p. 111732, https://doi.org/10.1016/j.jcp.2022.111732.
  * [53] N. Yang, H. Yang, and C. Yang, Multilevel field-split preconditioners with domain decomposition for steady and unsteady flow problems, Computer Physics Communications, 280 (2022), p. 108496, https://doi.org/10.1016/j.cpc.2022.108496.
  * [54] Y. Yang, S. Fu, and E. T. Chung, A two-grid preconditioner with an adaptive coarse space for flow simulations in highly heterogeneous media, Journal of Computational Physics, 391 (2019), pp. 1–13, https://doi.org/10.1016/j.jcp.2019.03.038.