On adaptive BDDC for the flow in heterogeneous porous media

TL;DR

This paper develops an adaptive BDDC method for simulating single-phase flow in heterogeneous porous media, improving computational efficiency through adaptive constraint selection and demonstrating effectiveness on benchmark problems.

Contribution

The paper introduces an adaptive algorithm for selecting flux constraints within BDDC, enhancing solver efficiency for flow in heterogeneous porous media.

Findings

Adaptive BDDC reduces iteration counts significantly.

Numerical upscaling properties observed in the first two steps.

Method performs well on 2D and 3D benchmark problems.

Abstract

We study a method based on Balancing Domain Decomposition by Constraints (BDDC) for a numerical solution of a single-phase flow in heterogenous porous media. The method solves for both flux and pressure variables. The fluxes are resolved in three steps: the coarse solve is followed by subdomain solves and last we look for a divergence-free flux correction and pressures using conjugate gradients with the BDDC preconditioner. Our main contribution is an application of the adaptive algorithm for selection of flux constraints. Performance of the method is illustrated on the benchmark problem from the 10th SPE Comparative Solution Project (SPE 10). Numerical experiments in both 2D and 3D demonstrate that the first two steps of the method exhibit some numerical upscaling properties, and the adaptive preconditioner in the last step allows a significant decrease in number of iterations of conjugate gradients at a small additional cost.

Tables

Table 1. Table 1: Substructuring of the 2D and 3D problems: N 𝑁 N is the number of subdomains, n Γ subscript 𝑛 Γ n_{\Gamma} is the number of (flux) degrees of freedom on interfaces, n f subscript 𝑛 𝑓 n_{f} is the number of faces, and n c subscript 𝑛 𝑐 n_{c} is the number of (initial) coarse degrees of freedom.

type of partitioning	$N$	$n_{\Gamma }$	$n_f$	$n_c$
2D
regular ($H/h=30$)	14	580	19	33
regular ($H/h=10$)	132	2360	236	368
irregular A	16	756	29	70
irregular B	64	1746	152	315
3D
regular ($H/h=10$)	27	5400	54	81
irregular	32	7267	129	335

Table 2. Table 2: Convergence of the non-adaptive method for the homogeneous case ( k = 1 𝑘 1 k=1 ) and the six layers of the SPE 10 problem. Here ϵ 0 subscript italic-ϵ 0 \epsilon_{0} and ϵ ∗ superscript italic-ϵ ∗ \epsilon^{\ast} are the errors (in % percent \% ) defined by ( 40 ), ω ~ ~ 𝜔 \widetilde{\omega} is the condition number indicator from ( 36 ), i ​ t 𝑖 𝑡 it  is the number of iterations for relative residual tolerance 10 − 6 superscript 10 6 10^{-6} , and κ 𝜅 \kappa  is the approximate condition number computed from the Lánczos sequence in conjugate gradients.

layer	$H/h=30$		$H/h=10$		irregular A		irregular B
	$it$	$\kappa$	$it$	$\kappa$	$it$	$\kappa$	$it$	$\kappa$
($k=1$)	11	2.790	14	3.980	12	3.151	14	4.050
1	15	8.879	22	9.491	17	6.714	19	11.197
20	14	5.749	19	6.926	15	6.524	18	6.429
60	162	4564.1	513	26,359.3	244	11,272.6	292	7301.7
85	183	9310.7	446	24,492.8	208	7170.4	392	58,931.7

Table 3. Table 3: Convergence of the adaptive method for layer  1 1 1 of the SPE 10 problem with the irregular partitioning A. Here τ 𝜏 \tau is the target condition number from Algorithm  7 , ϵ 0 subscript italic-ϵ 0 \epsilon_{0}  and  ϵ ∗ superscript italic-ϵ ∗ \epsilon^{\ast} are the errors (in  % percent \% ) defined by ( 40 ), ω ~ ~ 𝜔 \widetilde{\omega} is the condition number indicator from ( 36 ), i ​ t 𝑖 𝑡 it  is the number of iterations for relative residual tolerance 10 − 6 superscript 10 6 10^{-6} , and κ 𝜅 \kappa  is the approximate condition number computed from the Lánczos sequence in conjugate gradients. With τ = ∞ 𝜏 \tau=\infty no adaptive constraints were used, and (ms) indicates use of the multiscale constraints.

$\tau$	${ϵ}_0[\%]$	${ϵ}^{\ast }[\%]$	$\tilde{\omega }$	$n_c$	$it$	$\kappa$
$\infty$	73.21	30.55	13.586	70	17	6.714
(ms)	72.55	27.15	-na-	121	15	5.998
10	71.86	29.11	8.404	73	16	6.231
5	70.77	18.89	4.765	81	13	5.517
3	69.19	11.64	2.997	104	11	2.842
2	69.23	9.42	1.970	153	8	1.915

Table 4. Table 4: Convergence of the adaptive method for layer  85 85 85 with the irregular partitioning B.

$\tau$	${ϵ}_0[\%]$	${ϵ}^{\ast }[\%]$	$\tilde{\omega }$	$n_c$	$it$	$\kappa$
$\infty$	69.34	42.11	59,491.702	315	392	58,931.700
(ms)	68.76	39.00	-na-	494	297	8931.930
10,000	69.34	42.11	9275.614	316	347	9170.830
1000	68.03	40.29	898.754	360	152	836.227
100	67.21	38.38	98.117	430	54	95.439
10	66.16	35.68	9.885	489	19	9.672
5	63.07	31.85	4.888	536	13	4.836
3	56.19	18.87	2.988	614	10	3.010
2	53.44	14.87	1.997	743	7	1.879

Table 5. Table 5: Convergence of the non-adaptive method for the homogeneous case ( k = 1 𝑘 1 k=1 ) and the two 3D cutouts of the SPE 10 problem from Figure  6 . The headings are same as in Table  2 .

layer	$H/h=10$		irregular part.
	$it$	$\kappa$	$it$	$\kappa$
$\left(k=1\right)$	25	17.099	35	22.029
1–30	779	49,075.600	1968	$1.096\times {10}^6$
56–85	3762	$2.576\times {10}^6$	5277	$3.676\times {10}^6$

Table 6. Table 6: Convergence of the adaptive method for layers 1–30 of the SPE 10 with irregular partitioning.

$\tau$	${ϵ}_0[\%]$	${ϵ}^{\ast }[\%]$	$\tilde{\omega }$	$n_c$	$it$	$\kappa$
$\infty$	98.41	86.05	$1.191\times {10}^6$	335	1968	$1.096\times {10}^6$
(ms)	98.29	86.05	-na-	571	1943	$1.079\times {10}^6$
100,000	98.39	85.95	94,328.862	349	1280	92,307.000
10,000	98.14	84.27	9862.559	514	514	10,512.200
1000	97.41	82.73	995.230	989	175	1014.150
100	92.93	72.63	97.989	1331	60	108.673
10	87.47	66.66	9.993	1617	18	11.711
5	85.82	65.46	4.985	1898	13	6.069
3	82.90	63.22	2.997	2331	9	3.007
2	81.62	62.81	2.000	2997	6	1.930

Table 7. Table 7: Convergence of the adaptive method for layers 56–85 of the SPE 10 with the regular partitioning.

$\tau$	${ϵ}_0[\%]$	${ϵ}^{\ast }[\%]$	$\tilde{\omega }$	$n_c$	$it$	$\kappa$
$\infty$	99.31	66.24	$3.435\times {10}^6$	81	3762	$2.576\times {10}^6$
(ms)	99.30	66.24	-na-	99	3735	$2.566\times {10}^6$
100,000	99.33	65.78	95,129.959	122	1267	93,040.500
10,000	98.63	62.19	9834.429	188	498	9487.840
1000	98.25	63.81	990.920	373	183	1200.070
100	97.26	64.08	99.793	766	59	124.896
10	91.79	49.39	9.965	1154	17	9.506
5	88.88	43.77	4.991	1342	13	5.545
3	86.52	41.75	2.997	1610	9	3.205
2	85.18	41.05	2.000	1960	6	1.918

Equations

$$k^{-1}\mathbf{u}+\nabla p$$

$$\nabla \cdot \mathbf{u}$$

$$p$$

$$\mathbf{u}\cdot \mathbf{n}$$

$$-{\int }_{\Omega }{f}_{\Omega }dx+{\int }_{\partial \Omega }{g}_{E}ds=0,$$

$${\mathbf{H}}_0(\Omega ;\mathrm{div})=\left\{\mathbf{v}:\mathbf{v}\in L^2(\Omega );\nabla \cdot \mathbf{v}\in L^2(\Omega )\text{and}\mathbf{v}\cdot \mathbf{n}=0\text{on }\partial \Omega \right\},$$

$${\Vert \mathbf{v}\Vert }_{{\mathbf{H}}_0(\Omega ;\mathrm{div})}^2={\Vert \mathbf{v}\Vert }_{L^2(\Omega )}^2+H_{\Omega }^2{\Vert \nabla \cdot \mathbf{v}\Vert }_{L^2(\Omega )}^2,$$

$$L_0^2\left(\Omega \right)=\left\{q:q\in L^2\left(\Omega \right)\text{and}{\int }_{\Omega }qdx=0\right\}.$$

$${\int }_{\Omega }{k}^{-1}\mathbf{u}\cdot \mathbf{v}dx-{\int }_{\Omega }p\left(\nabla \cdot \mathbf{v}\right)dx$$

$$-{\int }_{\Omega }\left(\nabla \cdot \mathbf{u}\right)qdx$$

$$a\left(u,v\right)$$

$$b\left(u,q\right)$$

$$⟨f,q⟩$$

$$a\left(u,v\right)+b\left(v,p\right)$$

$$b\left(u,q\right)$$

$$Q=Q_0\oplus Q_I,Q_I=Q^1\times \cdots Q^N,$$

$${\int }_{\Omega }{q}_{0}dx=0,\forall q_0\in Q_0\text{and}{\int }_{{\Omega }^i}{q}^{i}dx=0,\forall q^i\in Q^i,i=1,\dots ,N.$$

$$W=W^1\times \cdots \times W^N.$$

$$a\left(u_I,v_I\right)+b\left(v_I,p_I\right)$$

$$b\left(u_I,q_I\right)$$

$$U=\hat{W}\oplus U_I,\hat{W}{\perp }_a{U}_I.$$

$$b\left(u_I,q_0\right)=-{\int }_{\Omega }\left(\nabla \cdot u_I\right)q_{0}dx=0.$$

$${\int }_{{\Gamma }^{ij}}\left(w^i-w^j\right)\cdot {\mathbf{n}}^{i}ds=0,w^i\in W^i,w^j\in W^j$$

$$\tilde{W}={\tilde{W}}_{\Delta }\oplus {\tilde{W}}_{\Pi }.$$

$$\tilde{W}=\left\{w\in \left(I-P\right)W:C_U\left(I-E\right)w=0\right\}.$$

$$Q^i=\left\{\left(q\in Q\right){|}_{{\Omega }^i}:C_Q^{i}q=0\right\},$$

$$a\left({\tilde{w}}_0,{\tilde{v}}_{\Pi }\right)+b\left({\tilde{v}}_{\Pi },p_0\right)$$

$$b\left({\tilde{w}}_0,q_0\right)$$

$$u_0=E{\tilde{w}}_0.$$

$$a\left(u_I,v_I\right)+b\left(v_I,p_I\right)$$

$$b\left(u_I,q_I\right)$$

$$u^{\ast }=u_0+u_I.$$

$$a\left(u_{\text{corr}},v\right)+b\left(v,p\right)$$

$$b\left(u_{\text{corr}},q\right)$$

$$u=u^{\ast }+u_{\text{corr}}.$$

$$b\left(u^{\ast },q\right)=⟨f,q⟩,\forall q\in Q.$$

Full text

1 Introduction

The Balancing Domain Decomposition by Constraints (BDDC), proposed independently by Cros [7], Dohrmann [9], and Fragakis and Papadrakakis [16], is one of the most popular methods of iterative substructuring. The method was developed as a preconditioner for the solution of systems of linear equation obtained by finite element discretizations of elliptic problems, and it has been originally derived as a primal counterpart of the Finite Element Tearing and Interconnecting - Dual, Primal (FETI-DP) method by Farhat et al. [14, 15]. Over the years the BDDC has been extended to other types of problems, for example to the nearly incompressible elasticity by Dohrmann [10], the Stokes problem by Li and Widlund [26], or advection-diffusion problems by Tu and Li [25, 46]. It is also relatively straightforward to extend the BDDC into multiple levels, as noted by Dohrmann [9]. The three-level methods were developed in two and three dimensions by Tu [43, 44], and Mandel et al. [29] extended the algorithm into a multilevel method within a more general multispace BDDC setting. Another class of problems, important in the context of this paper, is the flow in porous media based on mixed and mixed-hybrid finite element discretizations. Probably the first domain decomposition methods of this class were proposed by Glowinski and Wheeler [17]. Their Method II was preconditioned using BDD by Cowsar et al. [6], using BDDC by Tu [42], and Šístek et al. [36] extended this methodology to flow in porous media with combined mesh dimensions. This approach is regarded as hybrid because the method iterates on a system of dual variables (as Lagrange multipliers) enforcing the continuity of flux variables across the substructure interfaces. An alternative strategy, retaining the original primal variables was proposed by Tu [41, 45], who combined the BDDC preconditioner with an earlier algorithmic framework developed by Ewing and Wang [13], cf. also Mathew [31], which allows to solve the saddle-point problem obtained from mixed finite element discretization by conjugate gradients. The Nested BDDC by Sousedík [37] provided a multilevel extension by applying the framework from [41] recursively. Most recently, Zampini and Tu [49] presented another approach to multilevel BDDC including adaptive coarse space construction, which relies on a special, so-called, deluxe scaling.

There are two main ingredients of a BDDC method: a coarse space, which is defined by constraints on the values of degrees of freedom, and a scaling (averaging) operator, which provides a mapping between the solution space and the space in which the solves in the preconditioner are performed. The algorithm for adaptive selection of constraints for both methods the BDDC and FETI-DP was originally proposed by Mandel and Sousedík [28]. The algorithm was later generalized in a joint work with Šístek [30] into three spatial dimensions and implemented for the BDDC using an approach inspired by a partial subassembly and a change of variables by Li and Widlund [27]. Finally, we also reformulated the algorithm to treat the coarse space explicitly [38]. We note that there are many other approaches to the adaptive construction of the coarse spaces in BDDC, see [34] and the references therein, as well as for BDD see, e.g., [39]. There have been several scalings studied in the literature. In the multiplicity scaling, the weights are chosen proportionally to the number of subdomains sharing a given degree of freedom, and it is regarded as not robust for coefficient jumps. The $\rho$-*scaling *leads to robustness, but it relies on knowledge of the problem coefficients [22]. The stiffness scaling is based on the diagonal of the stiffness matrix, but in some cases with irregular meshes it may lead to high condition numbers [33, 35]. All these scalings involve diagonal matrices. Finally, the deluxe scaling introduced in [11] uses dense matrices, which are computed from inverses of localized Schur complements. It has been observed to be quite robust [32, 49] but also computationally intensive.

In this paper, we build on the primal strategy. The starting point is the two-level algorithm from [37], which we combine with adaptive selection of constraints following [38] and apply it to flow in heterogenous porous media. To this end, we use a reservoir from the 10th SPE Comparative Solution Project (SPE 10) cf., e.g., [1, 5] as the benchmark problem. The BDDC method from [37] solves for both flux and pressure variables. The fluxes are resolved in three steps: the coarse solve is followed by mutually independent subdomain solves, and last we look for a divergence-free flux correction and pressure using conjugate gradients (CG) with the BDDC preconditioner. The coarse solve in the first step is exactly the same as the coarse solve used in the BDDC preconditioner in the step three. It is assumed that the initial constraints preserve the iterates in a balanced subspace, in which the preconditioned operator is positive definite. Our goal here is to adapt the method to flow in realistic reservoirs, characterized by highly heterogeneous permeability coefficients in as simple way as possible. In particular, we translate the ideas used for elliptic problems in [38] to mixed formulations of flow in porous media discretized by the lowest-order Raviart-Thomas finite elements (RT0). The main component of the extension is the use of additional adaptive flux coarse basis functions. The starting point is the condition number bound formulated as a generalized eigenvalue problem, which is replaced by a number of local eigenvalue problems formulated for pairs of adjacent subdomains, and the eigenvectors, corresponding to the eigenvalues larger than a target condition number are used to construct the additional flux coarse basis functions. We note that from this perspective our method can be viewed as a way of numerical upscaling via the coarse basis functions known from the BDDC. Unlike [49] we do not use a change of basis and partial assembly of operators, and we also illustrate that for this problem the multiplicity scaling in combination with the adaptive algorithm and a simple diagonal rescaling of the pressure block in the setup of the problem is sufficient to construct a robust algorithm. Numerical experiments in both 2D and 3D demonstrate that the first two steps of the method exhibit some numerical upscaling properties, and the convergence rate of conjugate gradients in the last step can be estimated a priori in the setup of the adaptive algorithm.

The paper is organized as follows. In Section 2 we introduce the model problem, in Section 3 we recall the BDDC method and the preconditioner, in Section 4 we formulate the algorithm for adaptive selection of the flux constraints, in Section 5 we discuss some details of implementation, in Section 6 we present results of numerical experiments, and finallly in Section 7 we summarize and conclude our work.

For convenience, we identify finite element functions with the vectors of their coefficients in the corresponding finite element basis. These coefficients are also called variables or degrees of freedom. At a few places we will also identify linear operators with their matrices, in bases that will be clear from the context. For a symmetric positive definite bilinear form $a$, we will denote the energy norm by $\left\|u\right\|_{a}=\sqrt{a\left(u,u\right)}$.

2 Model problem

Let $\Omega$ be a bounded domain in $\mathbb{R}^{n}$, where $n=2$ or $3$. We would like to find the solution of the following mixed problem, which combines the Darcy’s law relating flux $\mathbf{u}$ and pressure $p$, and the equation of continuity,

[TABLE]

where $\partial\Omega=\overline{\Gamma}_{E}\cup\overline{\Gamma}_{N}$, and $\mathbf{n}$ denotes the unit outward normal of $\Omega$. The coefficient $k=k_{p}/\mu$, where $k_{p}$ is the permeability of the porous medium and $\mu$ is the viscosity of the fluid. For simplicity, we will set $\mu=1$ and so$~{}k=k_{p}$. Without loss of generality we will also assume that $\Gamma_{N}=\emptyset$, which requires a compatibility condition

[TABLE]

and the pressure $p$ will be uniquely determined up to an additive constant. We will further assume that $g_{E}=0$. These assumptions motivate the definition of a space

[TABLE]

equipped with the norm

[TABLE]

where $H_{\Omega}$ is the characteristic size of $\Omega$, and the definition of a space

[TABLE]

The weak form of the problem we wish to solve, is

[TABLE]

We refer, e.g., to the monographs [4, 40] for additional details and discussion.

Next, let $U$ be the lowest-order Raviart-Thomas (RT0) finite element space with a zero normal component on $\partial\Omega$ and let $Q$ be a space of piecewise constant finite element basis functions with a zero mean on $\Omega$. These two spaces, defined on the triangulation$~{}\mathcal{T}_{h}$ of$~{}\Omega$, where $h$ denotes the mesh size, are finite dimensional subspaces of $\mathbf{H}_{0}(\Omega;\operatorname{div})$ and $L_{0}^{2}(\Omega)$, respectively, and they satisfy a uniform inf-sup condition, see [4]. Let us define the bilinear forms and the right-hand side by

[TABLE]

In the mixed finite element approximation of problem (6)–(7), we would like to find a pair of fluxes and pressures $\left(u,p\right)\in\left(U,Q\right)$ such that

[TABLE]

We note that $Q$ is a finite-dimensional subspace of $L_{0}^{2}\left(\Omega\right)$ and therefore the unique solvability of the mixed problem (11)–(12) is guaranteed.

In the next section, we will describe the components of the two-level Nested BDDC, which allows an efficient iterative solution of problem (11)–(12).

3 The BDDC method

Let us consider a decomposition of$~{}\Omega$ into a set of nonoverlapping subdomains $\Omega^{i}$, $i=1,\ldots,N,$ also called substructures, forming a quasi-uniform triangulation of $\Omega$ and denote the characteristic subdomain size by $H$. Each substructure is a union of finite elements with a matching discretization across the substructure interfaces. Let $\Gamma^{i}=\partial\Omega^{i}\backslash\partial\Omega$ be the set of boundary degrees of freedom of a substructure$~{}\Omega^{i}$ shared with another substructure$~{}\Omega^{j}$, $j\neq i$, and define the interface by $\Gamma=\cup_{i=1}^{N}\Gamma^{i}$. Let us define a face as an intersection $\Gamma^{ij}=\Gamma^{i}\cap\Gamma^{j}$, $i\neq j$ and let us denote by $\mathcal{F}$ the set of all faces between substructures. Note that with respect to the RT0 discretization we define only faces, but no corners (nor edges in 3D) known from other types of substructuring.

We will solve problems similar to (11)–(12) on each substructure. As we have noted, such problems determine the pressure uniquely up to a constant, so we consider the decomposition of the pressure space

[TABLE]

where $Q_{0}$ consists of functions that are constant in each subdomain and have a zero average over the whole domain$~{}\Omega$, and the product space$\ Q_{I}$ consists of functions that have zero weighted average over one subdomain at a time. That is,

[TABLE]

Next, let $W^{i}$ be the space of flux finite element functions on a substructure $\Omega^{i}$ such that all of their degrees of freedom on $\partial\Omega^{i}\cap\partial\Omega$ are zero, and let

[TABLE]

Hence $U\subset W$ can be viewed as the subspace of flux functions from $W$ such that $u\cdot\mathbf{n}$ is continuous across substructure interfaces. Define $U_{I}\subset U$ as the subspace of flux functions such that $u\cdot\mathbf{n}$ is zero on the interface $\Gamma$, i.e., the space of “interior” flux functions, and let us also define a mapping $P:w\in W\longmapsto u_{I}\in U_{I}$ such that

[TABLE]

Functions from $\left(I-P\right)W$ will be called Stokes harmonic, cf. [40, Section 9.4.2].

Let $\widehat{W}$ be the space of Stokes harmonic functions that are continuous across substructure interfaces, and such that

[TABLE]

We note that from the divergence theorem, for all $u_{I}\in U_{I}$ and $q_{0}\in Q_{0}$, we obtain

[TABLE]

The BDDC is a two-level method characterized by a selection of certain coarse degrees of freedom. In the present setting these will be flux averages over faces shared by a pair of substructures at a time and pressure averages over each substructure. Let us denote by $\widetilde{W}\subset\left(I-P\right)W$ the subspace of Stokes harmonic functions such that their flux coarse degrees of freedom on adjacent substructures coincide; for this reason we will use the terms coarse degrees freedom and constraints interchangeably. Specifically, we define a zero-net flux constraint for a face$~{}\Gamma^{ij}$ as

[TABLE]

where $\mathbf{n}^{i}$ denotes the unit outward normal of $\Omega^{i}$.

Assumption 1**.**

Initial flux constraints (16) are prescribed over all faces.

This set of initial constraints will be enriched by the adaptive method described in Section 4. Now, let us define $\widetilde{W}_{\Pi}\subset\widetilde{W}$ as the subspace of functions with values given by the flux coarse degrees of freedom between adjacent substructures, and such that they are Stokes harmonic, and let us also define $\widetilde{W}_{\Delta}\subset\widetilde{W}$ as the subspace of all function such that their flux coarse degrees of freedom vanish. The functions in$~{}\widetilde{W}_{\Pi}$ are uniquely determined by the values of their coarse degrees of freedom, and

[TABLE]

The next ingredient is the projection $E:\widetilde{W}\rightarrow\widehat{W}$ defined by taking a weighted average of corresponding degrees of freedom on substructure interfaces, cf. Remark 3.

In implementation, we define$~{}\widetilde{W}$ using a matrix$~{}C_{U}$, which is a block diagonal with blocks$~{}C_{U}^{i}$, $i=1,\dots,N$, and it is constructed exactly as matrix $C$ in [28, Section 2.3],

[TABLE]

The values $C_{U}v$ will be called local flux coarse degrees of freedom, and the space$~{}\widetilde{W}$ consists of all functions such that their flux coarse degrees of freedom on adjacent substructures have zero jumps. The decomposition of the space $Q_{I}$ given by (13) can be also managed by constraints. We remark that this is somewhat non-standard practice in substructuring, because the constraints are commonly related only to the degrees of freedom at the interfaces. So, we define a space$~{}Q^{i}$, for $i=1,\dots,N$, as

[TABLE]

where the matrices $C_{Q}^{i}$ are selected so that (14) is satisfied. In implementation, $C_{Q}^{i}$ is a row vector with entries given by volumes of finite elements in subdomain $i$. Now we have all ingredients to recall the two-level BDDC method [37, Algorithm 2].

Algorithm 2** (BDDC method).**

Find the solution $\left(u,p\right)\in\left(U,Q\right)$ of problem (11)–(12) by computing:

the coarse component $u_{0}\in\widehat{W}$: solving $\left(\widetilde{w}_{0},p_{0}\right)\in\left(\widetilde{W}_{\Pi},Q_{0}\right)$ from

[TABLE]

dropping $p_{0}$, and applying the projection

[TABLE] 2. 2.

the substructure components $\left(u_{I},p_{I}\right)\in\left(U_{I},Q_{I}\right)$ from

[TABLE]

dropping $p_{I}$, and adding the solutions as

[TABLE] 3. 3.

the correction and the pressure $\left(u_{\text{corr}},p\right)\in\left(U,Q\right)$ from

[TABLE]

Specifically, use the CG method with the BDDC preconditioner defined in Algorithm 4, using the same setup of the coarse problem as in (20)–(21).

Finally, the flux variables are obtained as

[TABLE]

Remark 3**.**

The difference between problems (11)–(12) and (23)–(24) is that the latter problem has a vanishing second component, and therefore the correction $u_{\text{corr}}$ is divergence-free by (24). Also, we note that the initial flux constraints constructed according to (16) do not allow scaling weights in the scaling operator$~{}E$ to vary along the interface in order for $u^{\ast}$ to satisfy

[TABLE]

Therefore, in our numerical experiments, we use the multiplicity scaling unless the coefficient jumps are aligned with subdomain interfaces, see also [37, Remark 2].

The application of the BDDC preconditioner for the computation of $u_{corr}$ using two- resp. three-level method was studied by Tu [41, 45]. In [37], we applied Algorithm 2 recursively. Here, we will introduce a specific construction of the space$~{}\widetilde{W}_{\Pi}$ but before doing so, let us discuss Step 3 of Algorithm 2 in more detail.

The first step in substructuring is typically reduction of the problem to interfaces. In particular, problem (23)–(24) is reduced to finding $\left(\widehat{w},p_{0}\right)\in\left(\widehat{W},Q_{0}\right)$ such that

[TABLE]

where $f^{\ast}\in\widehat{W}^{\prime}$ is the reduced right-hand side. In implementation, the interiors are eliminated by the static condensation, problem (25)–(26) is solved iteratively, and the interiors $\left(u_{I},p_{I}\right)\in\left(U_{I},Q_{I}\right)$ are recovered in the post-correction. The key observation is, cf. [40, Section 9.4.2], that if we define a balanced subspace

[TABLE]

problem (25)–(26) becomes equivalent to the positive definite problem

[TABLE]

This observation justifies use of the CG method preconditioned by the BDDC provided that an initial guess is balanced, e.g., zero, and the outputs of the preconditioner are also balanced. It also implies that the iterates are effectively performed with the flux unknowns, and the pressure components$~{}p_{0}$ are resolved in the coarse correction of the preconditioner. The precise formulation of the two-level BDDC preconditioner for saddle-point problems follows. It is the reduced variant of [37, Algorithm 3].

Algorithm 4** (BDDC preconditioner).**

Define the preconditioner $\left(r_{B},0\right)\in\left(\widehat{W}^{\prime},Q_{0}^{\prime}\right)\longmapsto\left(\widehat{w},p_{0}\right)\in\left(\widehat{W},Q_{0}\right)$ by computing:

the coarse correction $\left(w_{\Pi},p_{0}\right)\in\left(\widetilde{W}_{\Pi},Q_{0}\right)$ from

[TABLE] 2. 2.

the substructure correction $w_{\Delta}\in\widetilde{W}_{\Delta}$ from

[TABLE] 3. 3.

the sum and average of the two corrections

[TABLE]

In order to state the condition number bound, we also need to introduce a larger space of balanced functions $\widetilde{W}_{B}$ such that $\widehat{W}_{B}\subset\widetilde{W}_{B}$ defined as

[TABLE]

The space$~{}\widetilde{W}_{\Pi}$ is also balanced, i.e., $\widetilde{W}_{\Pi}\subset\widetilde{W}_{B}$ by (28). Then also the output of the preconditioner (29) satisfies $\widehat{w}\in\widehat{W}_{B}$, and we refer to [37, Lemma 3] for the proof.

Finally, we formulate the condition number bound. If we note that $E$ is a projection, it is the same as [37, Theorem 4] or [41, Theorem 6.1], cf. also [28, Theorem 3].

Theorem 5**.**

The condition number $\kappa$ of the BDDC preconditioner from Algorithm 4 satisfies

[TABLE]

The bound$~{}\omega$ in (30) inspires the adaptive selection of the flux constraints.

4 Adaptive selection of the flux constraints

The basic idea is same as in our previous work on adaptive BDDC for elliptic problems [28, 30, 38]. The bound $\omega$ in (30) is equal to the maximal eigenvalue $\lambda_{\max}$ of the generalized eigenvalue problem

[TABLE]

From the Courant-Fisher-Weyl minimax principle cf., e.g., [8, Theorem 5.2], the bound$~{}\omega$ can be decreased by adding constraints in the definition of the space$~{}\widetilde{W}_{B}$ as:

Lemma 6** ([30, 38]).**

The generalized eigenvalue problem (31) has eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots\lambda_{n}\geq 0$. Denote the corresponding eigenvectors by$~{}w_{\ell}$. Then, for any $k=1,\dots,n-1$, and any linear functionals $L_{\ell}$, $\ell=1,\dots,k$,

[TABLE]

with equality if

[TABLE]

Because solving the global eigenvalue problem (31) is computationally expensive, we replace it by a collection of much smaller problems defined for all pairs of adjacent substructures, where a pair of substructures is adjacent if they share a face. All quantities associated with a pair of adjacent substructures $\Omega^{i}$ and $\Omega^{j}$ will be denoted by a superscript$~{}^{ij}$. In particular, we define $W^{ij}=W^{i}\times W^{j}$, and the local space $\widetilde{W}^{ij}$ of Stokes harmonic functions that satisfy the initial constraints at the face$~{}\Gamma^{ij}$ by

[TABLE]

We note that the space$~{}\widetilde{W}_{B}^{ij}$ is balanced, which is an implication of Assumption 1.

In these settings (31) becomes a local problem to find $w\in\widetilde{W}_{B}^{ij}$ such that

[TABLE]

The bilinear form $a^{ij}$ is associated on$~{}\widetilde{W}_{B}^{ij}$ with the Schur complement$~{}S^{ij}$ defined with respect to the interfaces$~{}\Gamma^{i}$,$~{}\Gamma^{j}$, and is positive-definite, cf. [41, Lemma 3.1].

Now we can proceed in the same way as in [38]. Let us denote by$~{}\mathcal{C}$ the matrix corresponding to $C_{U}^{ij}\left(I-E^{ij}\right)$. The orthogonal projection onto null$~{}\mathcal{C}$ is given by

[TABLE]

and we implement the local generalized eigenvalue problems (33) as

[TABLE]

which can be either solved using a dense eigenvalue solver [28] or eventually, since

[TABLE]

a subspace iterations such as the LOBPCG method [23], which runs effectively in the factorspace, could be also used. From (34), we wish the constraints to satisfy

[TABLE]

That is, we would add into the matrix $C_{U}^{ij}$ the rows

[TABLE]

but because by [38, Proposition 1] each row can be split as c_{\ell}^{ij}=\left[\begin{array}[c]{cc}c_{\ell}^{i}&-c_{\ell}^{i}\end{array}\right] and either half of$~{}c_{\ell}^{ij}$ is used to augment the matrices$~{}C_{U}^{i}$ and$~{}C_{U}^{j}$, see (39). We note that, due to the discretization using RT0 elements, the added rows are readily available in the form used in substructuring. The adaptive BDDC algorithm follows.

Algorithm 7** (Adaptive BDDC).**

Find the smallest $k$ for every two adjacent substructures to guarantee that $\lambda_{k+1}\leq\tau$, where $\tau$ is a given tolerance threshold (the target condition number), and add the constraints (35) to the definition of $\widetilde{W}$.

After the adaptive constraints are added, we define the heuristic condition number indicator as the largest eigenvalue $\omega^{ij}$ of all local eigenvalue problems (33), that is

[TABLE]

Remark 8**.**

It has been shown in [49, Theorem 4.3], see also [34, Theorem 3.10] and [32, Theorem 3.3], that the condition number$~{}\kappa$ of the adaptive BDDC operator satisfies

[TABLE]

where $N_{F}$ is the maximum number of faces of any subdomain. We note that this bound is pessimistic due to the factor$~{}N_{F}^{2}$, and in fact we observed $\kappa\approx\widetilde{\omega}$ in all experiments.

5 Implementation remarks

First, we describe a rescaling used to preserve numerical stability of the method with highly heterogeneous permeability coefficients. The variational problem (11)–(12) can be written in the matrix form as

[TABLE]

Assuming that the mesh size$~{}h\approx 1$, the entries in$~{}A$ are$~{}O\left(k^{-1}\right)$ and the entries in$~{}B$ are$~{}O\left(1\right)$. In particular, in the case of the SPE 10 data set we get $k^{-1}\approx 10^{6}-10^{12}$, and we found that some of the subdomain matrices and the matrix of the coarse problem may appear numerically singular. Due to the discontinuous approximation of the pressure, $B$ is a block-diagonal rectangular matrix. Each block corresponds to a particular subdomain, and it can be rescaled, e.g., by an average of the diagonal entries of $A$ corresponding to the degrees of freedom in this subdomain. Collecting this scaling coefficients in a diagonal matrix$~{}D$, we replace (37) by

[TABLE]

and the pressure is recovered at the end of computations as $p=D\overline{p}$.

5.1 Coarse degrees of freedom

The selection of the flux coarse degrees of freedom or, equivalently, flux constraints entails construction of the matrix $C_{U}$ in the definition of the space $\widetilde{W}$ by (18). Similarly, the selection of the pressure constraints, which facilitate the decomposition (13), entails construction of the matrices $C_{Q}^{i}$, $i=1,\dots,N$, in the definition of the spaces$~{}Q^{i}$ by (19). Following the standard practice in substructuring, in implementation we work with global and local degrees of freedom and the corresponding spaces, and vectors from these spaces are related by a restriction operator (a zero-one matrix). Therefore, the matrix$~{}C_{U}$ is constructed as a block-diagonal matrix using blocks$~{}C_{U}^{i}$ that select local flux coarse degrees of freedom from all degrees of freedom of substructure $i$, see [28, Section 2.3] for details. In the mixed finite element settings, each local coarse degrees of freedom selection matrix is constructed simply by augmenting the matrix$~{}C_{U}^{i}$ by a row$~{}C_{Q}^{i}$ as

[TABLE]

and the matrices$~{}C_{U}^{i}$ may be further augmented by the adaptive algorithm, see (35).

5.2 Solution of the local generalized eigenvalue problems

The choice of an eigensolver for the eigenvalue problems (34) is a delicate one. In general, the decision whether to use a dense or a sparse eigensolver depends on the type of the eigenvalue problem, size of the substructures, dimension of the problem, availability of a preconditioner for a sparse solver, and conditioning and numerical sensitivity of the underlying problem. All these factors will clearly affect the overall computational cost and performance of the method. We note that the formulation (34) allows to use a matrix-free iterative method such as the LOBPCG [23] in the same way as for elliptic problems, including that it can be further preconditioned by a local version of the BDDC as suggested in [38, Section 5], see also [21]. However, we found that dense eigenvalue solvers are more suitable for the SPE 10 dataset due to their robustness, and we used Matlab function eig in the numerical experiments.

5.3 Computational cost

Clearly, the two most computationally expensive parts of the method are the setup of the constraints by solving the set of the local eigenvalue problems, and the factorization of the coarse problem. There are many eigenvalue problems to be solved, but they are small and can be solved in parallel—this feature is similar to the setup of multiscale finite element methods [12]. Assuming that these can be solved efficiently, the bottleneck in computations is the factorization of the coarse problem. Specifically, it is crucial for the application of the method to appropriately balance the effort in the preconditioner and the global linear solver through a judicious choice of$~{}\tau$. This could be, for example, achieved as follows: one can partition the domain into subdomains balancing the sizes of subdomains and assuming certain size of the coarse problem (and ideally also taking into account the coefficient jumps and minimizing the size of interfaces), solve the set of local eigenvalue problems, and based on the eigenvalues determine the number of additional adaptive constraints (and hence the value of$~{}\widetilde{\omega}$) which minimize the work needed to factor the coarse problem and the work needed by preconditioned conjugate gradients, including the coarse problem back-substitutions, needed to reduce the error to desired accuracy based on the well-known error reduction formula of conjugate gradients see, e.g., [18, Theorem 10.2.6].

6 Numerical experiments

We implemented the method in Matlab and studied its convergence for problems with large variations in the permeability coefficients $k$. In all experiments we used relative residual tolerance $10^{-6}$ as the convergence criterion for the conjugate gradients. First, we run a test with jumps in$~{}k$ aligned with substructure interfaces, see Figure 1. For this problem we used stiffness scaling, which is in case of the lowest-order Raviart-Thomas (RT0) elements equivalent to the $\rho$-scaling. This also implies that the stiffness scaling works well for irregular meshes (unlike for nodal elements). The conjugate gradients with the BDDC preconditioner converged in $15$ steps and the approximate condition number computed from the Lánczos sequence in conjugate gradients was $\kappa=4.046$; with $k=1$ the method converged in $14$ steps and $\kappa=4.050$, see the rightmost column in Table 2. In the remaining experiments, we focused on problems with highly heterogeneous coefficients, and we used the multiplicity scaling. Specifically, we simulated flow in a porous media given by Model 2 of the 10th SPE Comparative Solution Project [5], which is publicly available on the Internet111http://www.sintef.no/Projectweb/GeoScale/Results/MsMFEM/SPE10/ and, in particular, we used a Matlab dataset described in [1]. The dimensions of the full model are $1200\times 2200\times 170$ (ft), and the distribution of the coefficients$~{}k$ is given over a regular Cartesian grid with $60\times 220\times 85$ grid-blocks. We used several layers and two 3D cutouts of the model for our numerical experiments. For the experiments in 2D, we used layers 1, 20, 60 and 85 shown in Figures 2–3. In the top layers 1 and 20 the permeability is relatively smooth, whereas the bottom layers 60, and 85 are fluvial and they are characterized by a spaghetti of narrow high-flow channels. In all layers the permeabilities range over at least six orders of magnitude. To drive a flow, we impose an injection (source) and a production well (sink) in the lower-left and upper-right corners, respectively. The discretization of each layer by the quadrilateral RT0 finite elements yields $39,880$ degrees of freedom. The layers were partitioned into subdomains in four ways: using two geometrically regular partitionings with the coarsening ratios $H/h=30$ and $H/h=10$, and two irregular partitionings. The details of the partitionings are summarized in Table 1 and illustrated by Figures 2–3. For the experiments in 3D, we used two domains consisting of $30\times 30\times 30$ elements extracted from layers $1$–$30$ and $56$–$85$ of the SPE 10 problem shown in Figure 6. To drive a flow, we impose an injection (source) and a production well (sink) in two distant corners of the domain. The discretization by the hexahedral lowest-order Raviart-Thomas (RT0) finite elements yields $110,700$ degrees of freedom. The domain was partitioned into subdomains in two ways: using one geometrically regular partitioning with the coarsening ratio $H/h=10$, and an irregular partitioning. The details of the partitionings are summarized in Table 1 and illustrated by Figure 7. All irregular partitionings were obtained using METIS 4.0 [20], and in order to test the adaptive algorithm we did not take into account the permeability coefficients.

It is interesting to note that the adaptive flux coarse basis functions capture to some extent features of the solution on the finite element mesh, and the quality of this approximation improves as the threshold$~{}\tau$ in Algorithm 7 decreases. We illustrate this fact by relative errors of solutions $u_{0}$ and $u^{\ast}$ obtained in Steps 1 and 2 of Algorithm 2 with respect to the exact solution$~{}u_{\text{exact}}$ obtained by a direct solve of the full problem. Specifically, the two relative errors are reported in tables as

[TABLE]

We also compare the adaptive method with constraints inspired by Multiscale mixed finite element method (MsMFEM) cf. [12, Algorithm 2.5.2] or [2, Section 3.2.1]. In particular, instead of the local eigenvalue problems we solved local Darcy’s flow problems, that is local counterparts of problem (1)–(2), with the source term

[TABLE]

and zero flux boundary condition on $\partial\Omega^{i}\cap\partial\Omega^{j}$. The source distribution function is set to $w_{i}\left(x\right)=1/\left|\Omega^{i}\right|$ in all subdomains except those containing a well, in which

[TABLE]

to ensure a conservative approximation on the fine grid. In the numerical experiments we then used the set of basic constraints (16) enriched by solving the above problem and taking the values of flux degrees of freedom on $\partial\Omega^{i}\cap\partial\Omega^{j}$ as additional constraints. Nevertheless, we note that there are other more advanced solvers based on multiscale strategies available in the literature see, e.g., Yang et al. [48] or la Cour Christensen et al. [24], and a thorough comparison of the methods would be of independent interest.

The results of numerical experiments in 2D are summarized in Tables 2–4. Table 2 shows performance of the nonadaptive method for a homogeneous case with$~{}k=1$ and the layers of the SPE 10 problem. It can be seen that for layers$~{}1$, and $20$ the convergence does not significantly depend on the partitioning and it is also quite comparable to the homogeneous case with $k=1$. On the other hand, for layers $60$ and $85$ the variations in coefficients aggravate convergence, which is also quite sensitive to the partitioning. This holds, in particular, for layer$~{}60$ which contains both regions that are highly heterogeneous and relatively homogeneous. It can be also seen by comparing left and right columns in Table 2 that increasing the number of subdomains (that is decreasing the coarsening ratio$~{}H/h$) leads to higher condition numbers and increase in iteration counts for both regular and irregular partitionings. This is not the case in the standard theory of domain decomposition methods, but here we suspect it can be attributed to the jumps in coefficients and larger interfaces. The performance of the adaptive algorithm is illustrated by Tables 3–4. Table 3 shows convergence for layer $1$ with irregular partitioning A, and Table 4 shows convergence for layer $85$ with irregular partitioning B. It can be seen that in all cases lower values of the threshold$~{}\tau$ lead to fewer iterations, and the value of the condition number indicator$~{}\widetilde{\omega}<\tau$ is in a good agreement with$~{}\kappa$, which is the approximate condition number estimate obtained from the Lánczos sequence in conjugate gradients. The adaptive constraints also lead to more significant improvement in convergence than the multiscale constraints. The problem for layer$~{}85$ is particularly interesting. From the right panel in Figure 3 we see that the coefficient jumps have very large variations even on the interfaces, which can be seen in the left panel of Figure 4. The right panel displays the eigenvalues of the corresponding eigenproblem: $\lambda_{1}\approx 3769.5$ and all other eigenvalues are less than$~{}20$.  Figure 5 then displays $300$ largest eigenvalues of the (global) BDDC preconditioned operator without adaptivity and with adaptive BDDC and target condition number $\tau=100$. We see that without adaptivity there is a single largest eigenvalue: specifically $\lambda_{1}=59,492$ and $\lambda_{2}=9,258$. For the adaptive BDDC with $\tau=100$ we get $\lambda_{1}=96.3$. Comparing this plot with Table 4 we see that the adaptive BDDC with $\tau=100$ introduces $115$ adaptive constraints, which corresponds to the number of the largest eigenvalues removed from the spectrum of the BDDC preconditioned operator. We also note that adding a single adaptive constraint reduces the iteration count from $392$ to $347$, which corresponds to the large gap in the spectrum of the operator without adaptivity. Setting $\tau$ to a lower value, for example, $\tau=3$, roughly doubles the number of constraints and the number of iterations is reduced to approximately $10$. Also, the lower value of$~{}\tau$ improve the approximation quality of the first two steps of Algorithm 2 and, for example, with $\tau<3$ we get the error $\epsilon^{\ast}<20\%$.

The results of numerical experiments in 3D are summarized in Tables 5–7. It can be seen from Table 5 that the numbers of iterations are significantly higher than in 2D, and the convergence is slower for the fluvial bottom layers $56$–$85$ comparing with the relatively smooth top layers $1$–$30$. The increase in iterations becomes even more pronounced in the case of the irregular partitioning also due to larger interfaces. The results of experiments with the adaptive algorithm are summarized in Tables 6–7. As in the 2D case, lower values of the threshold$~{}\tau$ lead in all cases to fewer iterations, and the values of$~{}\tau$, $\widetilde{\omega}$ and$~{}\kappa$ are in close agreement. Again, the multiscale constraints provide only a slight improvement of convergence. Table 6 shows convergence for layers$~{}1$–$30$. It can be seen that despite higher condition number of the problem corresponding to the irregular partitioning, the adaptive algorithm leads allows to decrease the iteration counts for lower values of$~{}\tau$. As in 2D, the first few adaptive constraints allow to decrease the number of iterations by a fairly large amount: here adding $14$ constraints reduces the number of iterations from the initial value $1968$ to $1280$. However, for example with$~{}\tau=10$ the number of iterations decreases to $18$, however the number of constraints grows rather significantly from $335$ to $1617$. Finally, the values of$~{}\epsilon_{0}$ and$~{}\epsilon^{\ast}$ are quite larger compared to the 2D experiments. Table 7 shows convergence for layers$~{}56$–$85$ and the regular partitioning $H/h=10$, and the trends are quite similar as in the previous case. That is, the adaptive algorithm allows to control the convergence of conjugate gradients, but the number of adaptive constraints is relatively high in particular for lower values of $\tau$. These trends are in agreement with the qualitative observations made from Figure 5.

7 Conclusion

We studied a method for solution of single-phase flow in heterogeneous porous media. We have, in particular, shown that the idea of adaptive BDDC, previously used for elliptic problems, can be also applied in the context of the BDDC method for mixed finite element discretizations using the lowest-order Raviart-Thomas finite elements, and that the adaptive method works well with the usual types of scaling used in substructuring. We illustrated that the resulting algorithm can be successfully applied for adaptive selection of the coarse flux degrees of freedom using several examples corresponding to the SPE 10 benchmark model. The effect of the adaptive construction of the flux coarse basis functions is twofold. First, the first two steps of the BDDC method provide some approximation properties with respect to the exact solution of the full problem, in particular in 2D. Second, the coarse problem provides a better preconditioner for conjugate gradients used in the third step. We also compared the adaptive constraints with constraints inspired by multiscale mixed finite element method, and we found that the adaptive constraints outperform the multiscale constraints.

Next, we experimented with different partitionings of the domains into substructures. While the adaptive method is able to overcome these issues in many cases, it is evident that a suitable partitioning makes the adaptive method more efficient. We note that development of optimal partitioning strategies is an open problem cf., e.g., [3, 47]. However, our experiments indicate that if it is not possible to find a suitable partitioning, the best strategy is to simply minimize the size of interfaces, which may be achieved by a simple geometric partitioning, see also [19].