An Observation On The Uniform Preconditioners For The Mixed Darcy Problem

TL;DR

This paper investigates uniform preconditioners for the mixed Darcy problem, demonstrating that pressure Schur complement preconditioners can be effective when hydraulic conductivity is small, using an operator preconditioning framework.

Contribution

It introduces a K- and h-uniform block preconditioner for the mixed Darcy problem based on operator preconditioning, addressing challenges in scaling and inf-sup conditions.

Findings

Pressure Schur complement preconditioners are effective for small hydraulic conductivity.

The proposed preconditioner is robust and uniform with respect to K and mesh size h.

Establishment of a K-uniform inf-sup condition is achieved.

Abstract

When solving a multi-physics problem one often decomposes a monolithic system into simpler, frequently single-physics, subproblems. A comprehensive solution strategy may commonly be attempted, then, by means of combining strategies devised for the constituent subproblems. When decomposing the monolithic problem, however, it may be that requiring a particular scaling for one subproblem enforces an undesired scaling on another. In this manuscript we consider the H(div)-based mixed formulation of the Darcy problem as a single-physics subproblem; the hydraulic conductivity, K, is considered intrinsic and not subject to any rescaling. Preconditioners for such porous media flow problems in mixed form are frequently based on H(div) preconditioners rather than the pressure Schur complement. We show that when the hydraulic conductivity, K, is small the pressure Schur complement can also be utilised for H(div)-based preconditioners. The proposed approach employs an operator preconditioning framework to establish a robust, K- and h-uniform block preconditioner. The mapping properties of the continuous operator are a key component in applying the theoretical framework point of view. As such, a main challenge addressed here is establishing a K-uniform inf-sup condition with respect to appropriately weighted Hilbert intersection- and sum-spaces.

Tables

Table 1. Table 1. Condition numbers of the operators ℬ 1 ​ 𝒜 subscript ℬ 1 𝒜 \mathcal{B}_{1}\mathcal{A} and ℬ 2 ​ 𝒜 subscript ℬ 2 𝒜 \mathcal{B}_{2}\mathcal{A} . (left) K 𝐾 K is constant. (right) K 𝐾 K is discontinuous given by ( 39 ). Condition number of using the exact Schur complement (i.e. not − Δ h subscript Δ ℎ -\Delta_{h} ) is shown in the braces.

		$2^{-2}$	$2^{-3}$	$2^{-4}$	$2^{-5}$
${ℬ}_1$	1	1.1	1.1	1.1	1.1
	$1^{-2}$	1.1	1.1	1.1	1.1
	$1^{-4}$	1.1	1.1	1.1	1.1
	$1^{-6}$	1.1	1.1	1.1	1.1
	$1^{-8}$	1.1	1.1	1.1	1.1
${ℬ}_2$	1	1.1(2.4)	1.1(1.4)	1.1(1.4)	1.1(1.4)
	${10}^{-2}$	2.3(2.3)	2.3(2.3)	2.3(2.3)	2.3(2.3)
	${10}^{-4}$	2.9(2.6)	3.1(2.6)	3.1(2.6)	3.0(2.6)
	${10}^{-6}$	3.0(2.6)	3.3(2.6)	3.5(2.6)	3.5(2.6)
	${10}^{-8}$	3.0(2.6)	3.3(2.6)	3.5(2.6)	3.6(2.6)

Table 2. Table 2. Condition numbers of the operators ℬ 1 ​ 𝒜 subscript ℬ 1 𝒜 \mathcal{B}_{1}\mathcal{A} and ℬ 2 ​ 𝒜 subscript ℬ 2 𝒜 \mathcal{B}_{2}\mathcal{A} with permeability matrix 𝐊 𝐊 {\bf K} in ( 39 ). Numbers in parentheses denote the condition number of ℬ 2 ​ 𝒜 subscript ℬ 2 𝒜 \mathcal{B}_{2}\mathcal{A} when using the exact Schur complement to approximation − ∇ ⋅ ⁡ ( 𝐊 ​ ∇ ) -\operatorname{\nabla\cdot}({\bf K}\nabla) .

	$\theta$		$2^{-3}$	$2^{-4}$	$2^{-5}$	$2^{-6}$
${ℬ}_1$	$0$	1	1.4	1.4	1.4	1.4
		${10}^{-2}$	1.4	1.4	1.4	1.4
		${10}^{-4}$	1.4	1.4	1.4	1.4
		${10}^{-6}$	1.4	1.4	1.4	1.4
		${10}^{-8}$	1.4	1.4	1.4	1.4
	$\pi /4$	1	1.4	1.4	1.4	1.4
		${10}^{-2}$	1.4	1.4	1.4	1.4
		${10}^{-4}$	1.4	1.4	1.4	1.4
		${10}^{-6}$	1.4	1.4	1.4	1.4
		${10}^{-8}$	1.4	1.4	1.4	1.4

Table 3. Table 3. Condition numbers of ( 40 ) with preconditioners ( 41 ) and ( 42 ). For ℬ 1 subscript ℬ 1 \mathcal{B}_{1} and ℬ 1 K subscript superscript ℬ 𝐾 1 \mathcal{B}^{K}_{1} only the results from direct solver (computing the entire spectrum) are included as iterative solver did not converge within 1500 iterations.

		$2^{-2}$	$2^{-3}$	$2^{-4}$
${ℬ}_1$	$1$	3.5	4.0	4.0
	${10}^{-1}$	14.3	18.2	19.0
	${10}^{-2}$	91	118	139
	${10}^{-3}$	396	731	$1.0\times {10}^3$
	${10}^{-4}$	739	$2.5\times {10}^3$	$5.2\times {10}^3$
	${10}^{-6}$	817	$3.6\times {10}^3$	$14\times {10}^3$
	${10}^{-8}$	818	$3.6\times {10}^3$	$15\times {10}^3$
${ℬ}_1^K$	$1$	3.5	4.0	4.0
	${10}^{-1}$	29.8	39.6	41.8
	${10}^{-2}$	559	744	893
	${10}^{-3}$	7$\times {10}^3$	$14\times {10}^3$	$20\times {10}^3$
	${10}^{-4}$	4$\times {10}^4$	$15\times {10}^4$	$32\times {10}^4$
	${10}^{-6}$	5$\times {10}^5$	$22\times {10}^5$	$89\times {10}^5$
	${10}^{-8}$	5$\times {10}^6$	$22\times {10}^6$	$90\times {10}^6$

Equations

$$\frac{1}{K}\mathbf{u}-\nabla p$$

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

$$\left[\begin{array}[]{cc}A&B^{T}\\ B&0\end{array}\right]\left[\begin{array}[]{c}x\\ y\end{array}\right]=\left[\begin{array}[]{c}c\\ d\end{array}\right],$$

$$\left[\begin{array}[]{cc}A^{-1}&0\\ 0&(BA^{-1}B^{T})^{-1}\end{array}\right]\quad\mbox{and}\quad\left[\begin{array}[]{cc}(A+B^{T}B)^{-1}&0\\ 0&X?\end{array}\right].$$

$$\left[\begin{array}[]{cc}\alpha A&B^{T}\\ B&0\end{array}\right]\left[\begin{array}[]{c}x\\ y\end{array}\right]=\left[\begin{array}[]{c}c\\ d\end{array}\right].$$

$$\mathcal{B}_{1}=\left[\begin{array}[]{cc}(\alpha A+\alpha B^{T}B)^{-1}&0\\ 0&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left(\frac{1}{\alpha}I\right)^{-1}}\end{array}\right].$$

$$\mathcal{B}_{2}=\left[\begin{array}[]{cc}(\alpha A+B^{T}B)^{-1}&0\\ 0&I{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}^{-1}}+(B(\alpha A)^{-1}B^{T})^{-1}\end{array}\right].$$

$$\mathcal{B}=\left[\begin{array}[]{cc}(K^{-1}-\nabla\operatorname{\nabla\cdot})^{-1}&0\\ 0&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}I^{-1}}+(-\operatorname{\nabla\cdot}(K\nabla))^{-1}\end{array}\right].$$

$${\Vert f\Vert }_{\alpha X}^2={\alpha }^2{\Vert f\Vert }_X^2.$$

$$X\cap Y=\left\{f|f\in X,\text{ and }f\in Y\right\},X+Y=\left\{f+g|f\in X,\text{ and }g\in Y\right\}.$$

$$\Vert z{\Vert }_{X\cap Y}^2=\Vert z{\Vert }_X^2+\Vert z{\Vert }_Y^2$$

$$\Vert z{\Vert }_{X+Y}^2=\underset{\begin{array}{c}z=x+y\\ x\in X,y\in Y\end{array}}{\inf }\Vert x{\Vert }_X^2+\Vert y{\Vert }_Y^2.$$

$${\left(X+Y\right)}^{\prime }=X^{\prime }\cap Y^{\prime }.$$

$$\Vert z{\Vert }_{X\cap Y}=\max \left(\Vert z{\Vert }_X,\Vert z{\Vert }_Y\right),\text{ and }\Vert z{\Vert }_{X+Y}=\underset{\begin{array}{c}z=x+y\\ x\in X,y\in Y\end{array}}{\inf }\Vert x{\Vert }_X+\Vert y{\Vert }_Y,$$

$$\frac{1}{2}\left({\Vert z\Vert }_X^2+{\Vert z\Vert }_Y^2\right)\le \max {\left(\Vert z{\Vert }_X,\Vert z{\Vert }_Y\right)}^2\le {\Vert z\Vert }_X^2+{\Vert z\Vert }_Y^2,$$

$$\underset{\begin{array}{c}z=x+y\\ x\in X,y\in Y\end{array}}{\inf }\Vert x{\Vert }_X^2+\Vert y{\Vert }_Y^2\le {\left(\underset{\begin{array}{c}z=x+y\\ x\in X,y\in Y\end{array}}{\inf }\Vert x{\Vert }_X+\Vert y{\Vert }_Y\right)}^2\le 2\underset{\begin{array}{c}z=x+y\\ x\in X,y\in Y\end{array}}{\inf }\Vert x{\Vert }_X^2+\Vert y{\Vert }_Y^2$$

$$T\in \mathcal{L}(X_1\cap X_2,Y_1\cap Y_2)\cap \mathcal{L}(X_1+X_2,Y_1+Y_2),$$

$${\Vert T\Vert }_{\mathcal{L}(X_1+X_2,Y_1+Y_2)}\le \max \left({\Vert T\Vert }_{\mathcal{L}(X_1,Y_1)},{\Vert T\Vert }_{\mathcal{L}(X_2,Y_2)}\right).$$

$$a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)$$

$$b(\mathbf{u},q)$$

$$a(\mathbf{u},\mathbf{v})=\left(K^{-1}\mathbf{u},\mathbf{v}\right)\text{mbox}andb(\mathbf{u},q)=\left(\nabla \cdot \mathbf{u},q\right).$$

$$\mathcal{A}\left[\begin{array}{c}\mathbf{u}\\ p\end{array}\right]:=\left[\begin{array}{cc}{K}^{-1}& -\nabla \\ \nabla \cdot & 0\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ p\end{array}\right]=\left[\begin{array}{c}\mathbf{f}\\ g\end{array}\right].$$

$$\mathbf{V}=K^{-1/2}{\mathbf{L}}^2\cap {\mathbf{H}}_0(\mathrm{div}),Q=L^2/\mathbb{R}+K^{1/2}(H^1/\mathbb{R}).$$

$$\mathbf{Z}=\{\mathbf{u}\in \mathbf{V}|b(\mathbf{u},q)=0,\forall q\in Q\}.$$

$$a(\mathbf{z},\mathbf{z})=\Vert \mathbf{z}{\Vert }_{K^{-1/2}{\mathbf{L}}^2}^2=\Vert \mathbf{z}{\Vert }_{\mathbf{V}}^2,\forall \mathbf{z}\in \mathbf{Z}$$

$$a(\mathbf{u},\mathbf{v})=(\mathbf{u},\mathbf{v}{)}_{K^{-1/2}}\le \Vert \mathbf{u}{\Vert }_{\mathbf{V}}\Vert \mathbf{v}{\Vert }_{\mathbf{V}}\forall \mathbf{u},\mathbf{v}\in \mathbf{V}.$$

$$b(\mathbf{v},q)$$

$$\text{SS}\in \mathcal{L}(Q^{\prime },\mathbf{V}),$$

$${\Vert \text{SS}\Vert }_{\mathcal{L}(Q^{\prime },\mathbf{V})}\text{ is independent of }K,$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left(\operatorname{\nabla\cdot}\SS g,\psi\right)=\left\langle g,\psi\right\rangle,\text{ for any }g\in Q^{\prime}\text{ and }\psi\in Q.}$$

$$\underset{\mathbf{v}\in \mathbf{V}}{\sup }\frac{\left(\nabla \cdot \mathbf{v},p\right)}{{\Vert \mathbf{v}\Vert }_{\mathbf{V}}}\ge \underset{g\in Q^{\prime }}{\sup }\frac{\left(\nabla \cdot \text{SS}g,p\right)}{{\Vert \text{SS}g\Vert }_{\mathbf{V}}}=\underset{g\in Q^{\prime }}{\sup }\frac{⟨g,p⟩}{{\Vert \text{SS}g\Vert }_{\mathbf{V}}}\ge C_S^{-1}\underset{g\in Q^{\prime }}{\sup }\frac{⟨g,p⟩}{{\Vert g\Vert }_{Q^{\prime }}}=C_S^{-1}{\Vert p\Vert }_Q.$$

$$\left(\nabla \varphi ,\nabla \psi \right)=-⟨g,\psi ⟩,\forall \psi \in H^1/\mathbb{R}.$$

Full text

Keywords. Uniform block preconditioner, Uniform inf-sup condition, Operator preconditioning framework, Mixed Darcy problem, Porous media, Hydraulic conductivity, Rescaling.

1. Introduction

In this paper we will consider the mixed formulation of the Darcy problem of the form

[TABLE]

equipped with suitable boundary conditions. The variables $\mathbf{u}$ and $p$ represent the fluid flux and pressure, respectively, and $K=\kappa/\mu_{f}$ denotes the hydraulic conductivity where $\kappa$ is the material permeability and $\mu_{f}$ the fluid viscosity.

In this manuscript we pursue two purposes. The first purpose is to establish a uniform-in-$K$ inf-sup condition for (1)-(2). Our second, and primary, purpose is the construction of some efficient block-diagonal preconditioners, for (1)-(2), exhibiting robustness for $K\in(0,1)$. These two objectives are connected. Indeed the operator preconditioning framework [15] uses the former stability result to provide for the latter a robust, block-diagonal preconditioner realization.

In general, the framework approach is predicated on establishing a well-posedness result for the continuous problem in $K$-weighted Sobolev spaces. For basic linear systems of the form

[TABLE]

there are, generally speaking, two common approaches for constructing block-diagonal preconditioners. One may utilize a Schur complement for the first unknown or, alternatively, the second unknown. That is, the structural options for the preconditioner are:

[TABLE]

The first approach results in three distinct, unit-sized eigenvalues  [16] for any $0<K<\infty$. On the other hand, to the authors knowledge, only partial explanations have been offered for the second approach for the mixed Darcy problem. In [1], ${\mathbf{H}}(\operatorname{div})$ preconditioners were constructed and applied to (1)-(2) in the case $K=1$ for which $X$ coincides with the inverse of a potentially diagonalized mass matrix. The same authors also developed multilevel methods for weighted ${\mathbf{H}}(\operatorname{div})$ spaces in [2], but did not discuss the corresponding appropriate scaling for a Darcy problem. In [17] the consequences of $K$-scaling, also of spatially varying $K$, were studied for both approaches; it was shown that the eigenvalues of the preconditioned system were affected significantly by $K$ and good results were obtained only when a proper rescaling was used.

If we now consider a scaled version (3) we arrive at a system with the general form

[TABLE]

In [11, 21, 22] it was suggested that problems of the form of (5) may be preconditioned efficiently by a scaling of the second structural option of (4); that is, preconditioners structurally based on

[TABLE]

This view was advanced in the aforementioned work, to develop preconditioners for the Biot and Brinkman problems, where the Darcy part was preconditioned by the approach of (6). The Brezzi conditions for the Darcy problem in the weighted norms corresponding to the choice of (6), that is $\frac{1}{\sqrt{K}}{\mathbf{H}}_{0}(\operatorname{div})\times\sqrt{K}L_{0}^{2}$, follow directly from the unscaled version by a simple scaling. Furthermore, the approach is closely related to the augmented Lagrangian approach which was investigated in [3, 9] and used successfully for the Oseen and Maxwell type problems, respectively.

Given the form of (6) we are motivated to ask the following question: is the $\alpha$-scaling of $B^{T}B$, in the top-left block, necessary? The motivation for this question comes from considering (1)-(2) as part of a multi-physics problem. For instance, as a single-physics problem, the isolated case of $K\rightarrow 0$, in (1)-(2), is not necessarily intrinsically interesting—a simple scaling resolves the issue of a vanishing $K$. However, in a multi-physics setting, solution algorithms are often constructed via decomposition into single-physics subproblems. Thus, requiring a certain scaling of one of the single-physics sub-problems, as we see here with the Darcy problem, may enforce an undesired scaling on other single-physics problems within the multi-physics system.

To investigate this question we consider an alternative to (6); namely

[TABLE]

In the current work we construct block preconditioners based on the operator preconditioning framework [1, 15]. In the continuous case, for the Darcy problem, the preconditioner takes the particular form, of (7), given by

[TABLE]

It will be shown that preconditioners based on $\mathcal{B}$ are robust with respect to both the mesh size and the permeability $K$. We will also consider permeabilites with jumps and anisotropy in our numerical experiments where we compare our proposed approach, (7), to the previous approach of (6). The remainder of this manuscript is organized as follows: Section 2 introduces the necessary notation and basic results; Section 3 discusses the continuous uniform stability (inf-sup) condition and the resulting continuous preconditioner; Section 4 addresses the corresponding discrete case of each; in Section 5 the preconditioners proposed in this paper are validated through numerical experiments. We also use this section to explore some cases of practical interest which are not covered directly by the theory of Sections 3 and 4. Finally, Section 6 offers some concluding remarks. In closing, we mention that the theory of Sections 3–4 assume $K\in(0,1)$ is an arbitrary, but fixed, constant.

2. Preliminaries

Let $\Omega$ be a bounded, connected Lipschitz domain in $\mathbb{R}^{n},\,n=2,3$. Then $L^{2}=L^{2}(\Omega)$ denotes the space of square integrable functions, whereas $H^{k}=H^{k}(\Omega)$ denotes the Sobolev space of functions with all derivatives up to order $k$ in $L^{2}$. The spaces $L^{2}/\mathbb{R}$ and $H^{1}/\mathbb{R}$ contain functions in $L^{2}$ and $H^{1}$, respectively, with zero mean value. Vector valued functions, and Sobolev spaces of vector valued functions, are denoted by boldface. The space ${\mathbf{H}}(\operatorname{div})$ contains functions in ${\bf L}^{2}$ with divergence in $L^{2}$ and the subspace of functions $\mathbf{u}\in{\mathbf{H}}_{0}(\operatorname{div})$ are those with zero normal trace. The notation $(\cdot,\cdot)$ is used for the $L^{2}$ inner product and analogously for vector fields. The norm corresponding to the $L^{2}$ inner product is expressed with the canonical double-bar $\left\|\cdot\right\|_{L^{2}}$. The duality pairing between a real Hilbert space $X$ and its dual $X^{\prime}$ is $\left\langle\cdot,\cdot\right\rangle$. Suppose that $X$ is a Hilbert space and that $\alpha>0$ is a fixed real value. Then we denote by $\alpha X$ the Hilbert space whose elements coincide with the elements of $X$ and with norm

[TABLE]

For two Hilbert spaces, $X$ and $Y$, we denote by $\mathcal{L}(X,Y)$ the space of bounded linear maps from $X$ to $Y$. In the subsequent analysis we employ both the intersection and sum of two Hilbert spaces $X$ and $Y$. These composite spaces are formally defined as follows: Let $X$ and $Y$ be two Hilbert spaces which are both subspaces of some larger Hilbert space. The intersection and sum space are defined, respectively, as:

[TABLE]

In this manuscript we will be concerned with the case where $X$ and $Y$ are Hilbert spaces. In this case $X\cap Y$ and $X+Y$ are also Hilbert spaces with respect to the norms:

[TABLE]

and

[TABLE]

In addition [4, Theorem 2.7.1] if $X\cap Y$ is dense in both $X$ and $Y$, then

[TABLE]

Remark 1**.**

In [4] $X$, $Y$, $X\cap Y$ and $X+Y$ are Banach spaces. Again, see [4, Lemma 2.3.1]. The norms of $X\cap Y$ and $X+Y$ are given explicitly by the alternative definitions

[TABLE]

which is different from the norms above (8) and (9). In our context, where we aim to derive preconditioners, the norms (8) and (9) are, however, more convenient. We therefore detail the equivalence between the two different definitions. In particular,

[TABLE]

and further

[TABLE]

Hence, the norms employed here, (8) and (9), are equivalent with the norms in [4]. Finally, we remark that in the more general Banach space setting (10) is understood as an identification under an isometry, c.f. [4, Theorem 2.7.1], with respect to the norms of (11). In the Hilbert case the map providing the identification (10) is, instead, an isomorphism with respect to the (equivalent) norms (8) and (9).

Now suppose that $\{X_{1},X_{2}\}$ and $\{Y_{1},Y_{2}\}$ are pairs of Hilbert spaces such that both elements of each pairing are subspaces of some larger Hilbert space. If $T$ is a bounded linear operator from $X_{i}$ to $Y_{i}$ for $i=1,\,2$, then

[TABLE]

and in particular

[TABLE]

3. Continuous stability and preconditioning

The weak formulation of (1)–(2) reads: Find $\mathbf{u}\in\mathbf{V},p\in Q$ such that

[TABLE]

where

[TABLE]

We recall that $K\in(0,1)$ is considered an arbitrary but fixed constant. The corresponding coefficient matrix reads,

[TABLE]

The primary result of this section can now be stated:

Proposition 1** (A uniform-in-$K$ continuous stability result).**

Consider the problem (13)-(14) where $\mathbf{u}$ and $p$ are chosen, respectively, in the spaces:

[TABLE]

Then the Brezzi conditions are satisfied uniformly in $K$.

Proof.

We begin with two comments: first, for applications to preconditioning, we are particularly interested in the existence of an inf-sup condition for $b(\mathbf{v},q)$ that is independent of the choice of $K\in(0,1)$. We will prove that such a condition holds for the pair of spaces $\mathbf{V}$ and $Q$. Second, we remark that the arguments below are similar to [14], although the Sobolev spaces and bilinear forms involved are different. We now verify all of the requisite Brezzi conditions for (13)–(14). Towards this end let

[TABLE]

Then clearly

[TABLE]

and hence coercivity of $a(\cdot,\cdot)$ over $\mathbf{Z}$ is established. Furthermore, the boundedness of $a(\cdot,\cdot)$ follows from

[TABLE]

The boundedness of $b$ follows from a decomposition argument. Let $q=q_{0}+q_{1}$, $q_{0}\in L^{2}/\mathbb{R}$ and $q_{1}\in H^{1}/\mathbb{R}$ then

[TABLE]

Taking the infimum over all decompositions of $q$ yields the desired bound.

To establish the uniform inf-sup condition we will demonstrate the existence of a linear operator $\SS$ which satisfies the following properties:

[TABLE]

Suppose first that such an operator $\SS$ exists satisfying conditions (17)–(19), and define the constant $C_{\SS}=\left\|\SS\right\|_{\mathcal{L}(Q^{\prime},\mathbf{V})}$. Under these assumptions the inf-sup condition follows directly. To see this, take any $p\in Q$; it follows that:

[TABLE]

We now show that such an operator $\SS$ can be realized as the solution of a suitable Poisson problem. Moreover, the $K$-independence of the operator norm will follow from a suitable scaling between the primal and mixed formulations of the problem. Towards this end let $g\in(H^{1}/\mathbb{R})^{\prime}$ and define $\phi\in H^{1}/\mathbb{R}$ to be the solution of the homogeneous Neumann problem, i.e.,

[TABLE]

Now define $\SS$ by $\SS g=\nabla\phi$; then (21) implies

[TABLE]

When interpreted in the weak sense, (21) provides the identity (19); that is $\operatorname{\nabla\cdot}\SS g=g$ where $\operatorname{\nabla\cdot}:\mathbf{L}^{2}\rightarrow(H^{1}/\mathbb{R})^{\prime}$.

For a second perspective on establishing the operator $\SS$, we now assume that $g\in(L^{2}/\mathbb{R})^{\prime}$. Consider then a mixed formulation of (13)–(14), with $K=1$, given by: find $(\mathbf{r},\phi)\in{\mathbf{H}}_{0}(\operatorname{div})\times L^{2}/\mathbb{R}$ such that

[TABLE]

The above problem is classical, and as a result of its well-posedness we define $\SS g=\mathbf{r}$. It is a straightforward exercise to establish that, for $g\in L^{2}/\mathbb{R}$, the above definition of $\SS g$ coincides with that of (21). Moreover we have that $\operatorname{\nabla\cdot}\SS g=g$, from the second equation of (23), and

[TABLE]

The primary observation, then, is that a scaling argument now reveals that $\SS\in\mathcal{L}(Q^{\prime},\mathbf{V})$ with operator norm independent of $K$. To see this, first scale (22) by $K^{-1/2}$; i.e.  (22) implies

[TABLE]

Combine the above with (24) and use (12) to get

[TABLE]

for any $K>0$.

In particular note that the operator norm of $\SS$ is independent of $K$. This follows since the operator norms related to (22), both before and after rescaling, and (24) are independent of $K$. We have shown that $\SS$, defined by (22)–(24), satisfies the properties (17)–(19). Thus, the desired uniform inf-sup condition is established, pursuant to the discussion preceding (20). ∎

Remark 2**.**

Note that in the proof of Proposition 2, above, we have made use of (10), and reflexivity, to establish that

[TABLE]

With a $K$-uniform stability established, the framework put forth in [15] details the construction of a robust and efficient preconditioner for $\mathcal{A}$ can be constructed from preconditioners for the operators realizing the $\mathbf{V}$- and $Q$-norms. The following proposition states precisely this:

Corollary 1** (A $K$-robust preconditioner for the continuous mixed Darcy problem).**

The preconditioner defined by

[TABLE]

provides a robust preconditioner for (13)–(14) in the sense that the condition number of $\mathcal{B}\mathcal{A}$ is bounded uniformly in $K$. Here, $\mathcal{A}$ is given by (15).

Proof.

Preconditioners for the $\mathbf{V}$ inner product are well-known, cf. e.g. [2, 7, 10, 12, 13]. Following the framework of [15] it therefore suffices to construct an operator realizing the $Q$-norm. The result can then be combined with the operator realizing the $\mathbf{V}$-norm from the sources cited above. For the $Q$-norm, we begin by remarking that (9) is equivalent to the expression:

[TABLE]

In fact, the infimum on the right hand side is attained by the unique $\phi\in H^{1}/\mathbb{R}$ solving the variational problem:

[TABLE]

The above corresponds to the statement that $\phi$ is the solution of the elliptic problem

[TABLE]

where $\mathbf{n}$ is the outward pointing normal to $\partial\Omega$, and recall that $q\in L^{2}/\mathbb{R}$. We will make use of the expressions $q-\phi=-K\Delta\phi$ and, equivalently, $\phi=(I-K\Delta)^{-1}q$. Since $\phi$ minimizes (26), we have from the observations above that

[TABLE]

The canonical preconditioner, $\mathcal{B}_{Q}:Q^{\prime}\to Q$, corresponding to the above is

[TABLE]

According to [15], a continuous and robust-in-$K$ preconditioner for (15), is then

[TABLE]

This is precisely (25).

∎

4. Discrete Stability and preconditioning

In this section we describe the construction of a preconditioner for discretizations based on Brezzi-Douglas-Marini (BDM) and Raviart-Thomas (RT) elements  [5, 18]. Again we assume $K\in(0,1)$ is an arbitrary but fixed constant. The discrete approach reflects many aspects of the continuous setting of section 3. However, due to the discontinuous polynomial nature of the pressure elements, we first define a discrete $H^{1}$-norm to establish the $Q$-norm in the discrete case.

Let $\mathcal{T}_{h}$ be a shape regular simplicial mesh defined on the bounded, Lipschitz domain $\Omega$ and let $r\geq 0$. Let $\mathbf{V}_{h}$ be the ${\mathbf{H}}(\operatorname{div})$-conforming discrete space given by either the RT elements of order $r$ or the BDM elements of order $r+1$. Define $Q_{h}$ to be the usual corresponding space of discontinuous, piecewise polynomials of order $r$. Consider the discrete mixed Darcy problem given by: find $\mathbf{u}_{h}\in\mathbf{V}_{h}$ and $p_{h}\in Q_{h}$ such that

[TABLE]

The discrete $Q_{h}$-norm will be defined in terms of a discrete gradient which is the negative $L^{2}$-adjoint of the $\operatorname{\nabla\cdot}$ operator on $\mathbf{V}_{h}$. First, $\nabla_{h}:Q_{h}\to\mathbf{V}_{h}$ is defined by

[TABLE]

It is well-known, [6, 8], that with these particular choices of $Q_{h}$ and $\mathbf{V}_{h}$, there is an $h$-independent constant $\beta>0$ such that

[TABLE]

for every $q\in Q_{h}$. It follows that $\nabla_{h}$ is injective and we can define the discrete $H^{1}$-norm on $Q_{h}$ via

[TABLE]

We denote the space $H_{h}^{1}$ as the set $Q_{h}$ equipped with the norm $\left\|\cdot\right\|_{1,h}$ and the space $L^{2}_{h}$ as the set $Q_{h}$ equipped with the usual $L^{2}$-norm. The discrete analogue of the $Q$-norm, i.e. the discrete form of (26) which is itself equivalent to the discrete form of (9) is given as

[TABLE]

The spaces ${\mathbf{H}}_{0,h}(\operatorname{div})$ and $\mathbf{L}^{2}_{h}$ are defined analogously to $L^{2}_{h}$ above. The primary result of this section, which we now state, is the discrete analogue of Proposition 2.

Proposition 2** (A uniform-in-$K$ discrete stability result).**

Consider the problem (29)-(30) where $\mathbf{u}_{h}$ and $p_{h}$ are chosen, respectively, in the spaces:

[TABLE]

Then, under the corresponding norms, the Brezzi conditions are satisfied uniformly in $K$.

Proof.

Boundedness and coercivity of $a(\cdot,\cdot)$ and boundedness of $b(\cdot,\cdot)$ follows from the same arguments put forth in section 3. Verifying a $K$-independent inf-sup condition will therefore conclude the argument. To accomplish this, a left-inverse of $\nabla_{h}$ will be constructed, satisfying appropriate bounds, allowing for a similar argument to that of section 3 for the operator $\SS$. Let $\mathbf{Z}_{h}$ denote the discrete kernel of the $\operatorname{\nabla\cdot}$ operator; i.e. the set of $\mathbf{v}_{h}\in\mathbf{V}_{h}$ for which

[TABLE]

From (32) it follows (cf. [8]) that $\nabla_{h}:Q_{h}\to\mathbf{Z}_{h}^{\perp}$ is a linear bijection. Furthermore, every $\mathbf{v}\in\mathbf{V}_{h}$ can be uniquely decomposed as

[TABLE]

where $r\in Q_{h}$, with $\nabla_{h}r\in\mathbf{Z}_{h}^{\perp}$, and $\tilde{\mathbf{v}}\in\mathbf{Z}_{h}$. Since the spaces considered for $\mathbf{V}_{h}$ satisfy the relation $\operatorname{\nabla\cdot}\mathbf{V}_{h}\subset Q_{h}$ it follows that $\operatorname{\nabla\cdot}\tilde{\mathbf{v}}=0$, for every $\tilde{\mathbf{v}}\in\mathbf{Z}_{h}$, and the decomposition (36) is orthogonal with respect to both the ${\mathbf{H}}(\operatorname{div})$- and $\mathbf{L}^{2}$ inner products.

We now define the lifting operator $\Theta_{h}:\mathbf{V}_{h}\to Q_{h}$ by $\Theta_{h}\mathbf{v}=r$, according to (36). It is evident that $\Theta_{h}\nabla_{h}$ is the identity operator on $Q_{h}$ and that $\Theta_{h}\tilde{\mathbf{v}}=0$ for all $\tilde{\mathbf{v}}\in\mathbf{Z}_{h}$. Moreover, the inf-sup condition (32) and the ${\mathbf{H}}(\operatorname{div})$-orthogonality of (36) implies that

[TABLE]

which means that $\Theta_{h}\in\mathcal{L}({\mathbf{H}}_{0,h}(\operatorname{div})^{\prime},L^{2}_{h})$. From the $L^{2}$-orthogonality of (36) we also have that $\left\|\Theta_{h}\mathbf{v}\right\|_{1,h}\leq\left\|\mathbf{v}\right\|$, which implies that $\Theta_{h}\in\mathcal{L}(\mathbf{L}_{h}^{2},H_{h}^{1})$. From these bounds on $\Theta_{h}$, together with (12), we deduce that

[TABLE]

Since $\Theta_{h}\nabla_{h}$ is the identity on $Q_{h}$, we get from (37) that for every $q\in Q_{h}$,

[TABLE]

where $C=\max(1,\beta^{-1})$ and is thus independent of both $h$ and $K$. ∎

Corollary 2** (A $K$ and $h$-robust preconditioner for the discrete mixed Darcy problem).**

Define $\mathcal{A}_{h}$ as the coefficient matrix characterizing the left-hand side of the discrete problem (29)–(30). Then the preconditioner defined by

[TABLE]

provides a robust preconditioner for (29)–(30). That is, the condition number of $\mathcal{B}_{h}\mathcal{A}_{h}$ is bounded uniformly with respect to both $K$ and $h$.

Proof.

The proof again follows the ideas of [15] and is nearly identical to the proof of Corollary 1. We therefore outline only the main ideas here. The arguments in the proof of Proposition 2 imply that $\mathcal{A}_{h}$ is a homeomorphism from $\mathbf{V}_{h}\times Q_{h}$ to its dual. Moreover, the norms on $\mathcal{A}_{h}$ and its inverse are bounded independently of $K$ and $h$. Arguments analogous to those defining $\mathcal{B}$, in Section 3, lead directly to the framework [15] preconditioner given by the operator $B_{h}:\mathbf{V}_{h}^{\prime}\times Q_{h}^{\prime}\to\mathbf{V}_{h}\times Q_{h}$ defined by (38).

∎

Remark 3**.**

We remark that a small liberty has been taken for the notation of $I_{h}$ in (38). Specifically, $I_{h}$ in the top left block of (38) signifies the identity $I_{h}:\mathbf{V}_{h}\rightarrow\mathbf{V}_{h}$ while the use of the same symbol in the bottom right block signifies the identity on $Q_{h}$; recall that the discrete gradient, $\nabla_{h}$, is defined by (31).

5. Numerical Experiments

Let $\Omega$ be a triangulation of the unit square such that the unit square is first divided in $N\times N$ squares of length $h=1/N$. Each square is then divided into two triangles. Below we will consider the case of homogeneous Dirichlet conditions for the flux; fluxes will be discretized by zero’th order RT elements. In addition, we approximate the pressure in the space of piecewise constants and compute the eigenvalues of the preconditioned system. In the sections that follow, three different paradigms are presented: in Section 5.1 a spatially constant permeability, pursuant to the theory of Section 3, is demonstrated; in Section 5.2 we consider a computational experiment extending a jump in the (scalar) coefficient, $K$, and an anisotropic conductivity tensor; finally, in Section 5.3, we compare the analogues of the preconditioning strategies analyzed here in the context of a coupled problem with Darcy subsystem.

5.1. A spatially constant conductivity

Order-optimal multilevel methods for both ${\mathbf{H}}(\operatorname{div})$ and $H^{1}$ problems are well-known; thus, we consider preconditioners based simply on exact inversion. It was shown in  [20] that the following local $-\Delta_{h}$ operator is spectrally equivalent to $-\operatorname{\nabla\cdot}\nabla_{h}$ for discontinuous Lagrange elements of arbitrary order:

[TABLE]

Here, $\mathcal{T}_{h}$ is a triangulation of the domain $\Omega$; internal faces are signified by the set $\mathcal{E_{I}}$ whereas $\mathcal{E_{D}}$ denotes faces at the boundary associated with a pressure Dirichlet condition. Furthermore $\left\{\!\left\{f\right\}\!\right\}=\tfrac{1}{2}(f|_{T^{+}}+f|_{T^{-}})$ and $[\![f]\!]=f|_{T^{+}}-f|_{T^{-}}$ are respectively the average and jump value of $f$ from the two elements $T^{\pm}$ that share the internal facet. In Table 1 we compare the discrete versions of the following two preconditioners

[TABLE]

and

[TABLE]

Table 1 shows that both $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ yield robust results for any $K\in(0,1)$, but that $\mathcal{B}_{1}$ is usually somewhat better. This is somewhat surprising since the opposite is the case in simple tests with random matrices. Numerical experiments (not reported here) confirm that the robust behavior applies also to Neumann conditions and BDM elements.

5.2. An anisotropic conductivity tensor with jump

We now consider the case where a jump is present in the (scalar) coefficient, $K$, and the conductivity tensor is given by an anisotropic matrix. These cases are not covered by the theoretical analysis and it is as such interesting to compare the alternative preconditioners. To this end let $\Omega$ be a unit square and

[TABLE]

where $K_{0}\in\left(0,1\right]$ and ${\bf R}_{\theta}$ is a matrix of rotation by angle $\theta$. Homogeneous Dirichlet conditions for the normal flux shall be imposed on the left and right edges. As in the previous experiments we consider a uniform triangulation of the domain and finite element discretization in terms of zero’th order RT and piecewise constant elements.

Varying the magnitude of the $K_{0}$ Table 1 shows that both preconditioners are practically unaffected by the presence of the considered discontinuity. We remark that the pressure preconditioner in the discrete $\mathcal{B}_{2}$ operator is discretized as follows

[TABLE]

Note in particular that $K$ is averaged using the harmonic mean. Table 2 shows the condition numbers of the preconditioned problems with matrix valued permeability. Here the preconditioners $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$ are generalized as

[TABLE]

with $k=1/\sum_{i}\lambda_{i}({\bf K})^{-1}$ and $\lambda_{i}({\bf K})$ the $i$-th eigenvalue of the permeability matrix.

Remark 4**.**

We note that the leading block of $\mathcal{B}_{1}$ is implemented as a solver for the variational problem. Namely, we find $\mathbf{u}\in{\mathbf{H}}_{0}(\operatorname{div})$ such that

[TABLE]

The preconditioner $\mathcal{B}_{1}$ is robust with respect to the conditioning of ${\bf K}$; in particular the simple scaling of the pressure mass matrix is sufficient. This stands in contrast to the case of preconditioner $\mathcal{B}_{2}$; here, the same scaling yields $h$ independent condition numbers only for relatively large $K_{0}$. In order to recover the mesh and parameter robustness attributes for $\mathcal{B}_{2}$ we approximate the operator $-\operatorname{\nabla\cdot}({\bf K}\nabla_{h})$ as the exact Schur complement of $\mathcal{A}_{h}$. That is, an exact inverse of ${\bf K}^{-1}I_{h}$ on $V_{h}$ is used in the construction of the preconditioner; c.f. Table 2. We remark that [17] discusses a construction based on a diagonalized mass matrix. In particular, for ${\bf K}$ diagonal, a Jacobi preconditioner is shown to yield bounds independent of ${\bf K}$.

5.3. A coupled problem with a Darcy subsystem

To illustrate the potential benefit of having the possibility to choose different scales of the pressure, we consider a simplified Biot problem. We remark that fully-parameter-robust preconditioners for Biot’s problem are established, in detail, in [11]; the discussion here is intended only to illustrate the difficulties encountered in a single parameter regime. With this caveat in mind: consider a simplified, three-field formulation of Biot’s problem that has been differentiated in time via an implicit Euler scheme; the resulting system has the form:

[TABLE]

where $\sigma(\mathbf{u})=2\epsilon(\mathbf{u})=\left(\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right)$ and $K\in(0,1]$ is assumed fixed, and spatially constant, but otherwise arbitrary. Here the system (40) is completed by the Dirichlet boundary conditions on the displacement $\mathbf{u}$, i.e. $\mathbf{u}=\mathbf{u}_{0}$ on $\Gamma_{D}\subset\partial\Omega$, and Neumann condition $\sigma\cdot\mathbf{n}=\mathbf{h}_{0}$ on $\Gamma_{N}=\partial\Omega\setminus\Gamma_{D}$. For the Darcy subproblem we then set $\mathbf{v}\cdot\mathbf{n}=\mathbf{v}_{0}$ on $\Gamma_{D}$ and $p=p_{0}$ on $\Gamma_{N}$. It is here clear that the pressure provides the coupling between the elastic deformation and the porous media flow and that Brezzi theory of the two sub-systems of Stokes and Darcy type enables stability of the coupled system [11, 19]. The Stokes problem requires that the displacement is bounded in $\mathbf{H}^{1}$, while the pressure is bounded in $L^{2}$. For the Darcy problem, a flux in ${\mathbf{H}}_{0}(\operatorname{div})$ can be readily combined with a pressure in $L^{2}$. However, when $K\rightarrow 0$, the norms must be scaled. An in particular, as proposed in [17, 21, 22] a natural norm for the pressure is $\sqrt{K}L^{2}$. However, the Stokes-type coupling between the displacement and pressure suggests the pressure should not be scaled with $K$. On the other hand, the pressure norm $L^{2}+K^{1/2}$, used in our paper, can be bounded by $L^{2}$. Below, we illustrate some of the problems that may occur in a preconditioning setting with wrongly scaled pressures and demonstrate that the scaling proposed in this paper can indeed be combined into a robust preconditioner for a simplified Biot problem.

As preconditioners for (40) we shall consider operators

[TABLE]

where, importantly, $\mathcal{B}_{2}$ avoids rescaling of the $\nabla\nabla\cdot$ block, cf. (6) and the related discussion of the issue in the context of multiphysics problems. In this respect, the operator is similar to the Darcy preconditioner (25). Our aim is then to demonstrate that such scaling may be not be desirable for $K$-robustness. To this end we perform a numerical experiment where $\Omega=(0,1)^{2}$ with $\Gamma_{D}$ the union of left, right and bottom edge of $\Omega$. Furthermore a uniform triangulation of the domain is considered with spaces of the discrete displacement, flux $\mathbf{v}$ and the pressure constructed respectively by continuous piecewise quadratic Lagrange elements, zero’th order RT elements and piecewise constant elements. Indeed, Table 3 shows that only the preconditioner $\mathcal{B}_{2}$ leads to condition numbers bounded in $K$. One may cautiously observe that the pressure block in (6) includes a $K$-scaling. It is a reasonable enquiry, then, to ask whether including this scaling can rectify the poor performance observed for the $\mathcal{B}_{1}$ preconditioner. To test this observation let us also consider preconditioning (40) with operators

[TABLE]

Condition numbers corresponding to these operators are shown in Table 3. It is evident that neither of the preconditioners of (42) present any improvement compared to their unscaled variants. We conclude, comparing $\mathcal{B}^{K}_{2}$ and $\mathcal{B}_{2}$, that rescaling the pressure is not necessary for robustness with respect to $K$.

6. Conclusion

In this paper we have introduced a new preconditioner for the Darcy problem and shown that for $K\in\mathbb{R}+$ the preconditioner is stable as $K\rightarrow 0$. The preconditioner is of the form

[TABLE]

We have compared it with a simpler approach where a scaling is introduced on the divergence term, i.e. the preconditioner of the form:

[TABLE]

The preconditioners were tested with both constant and small permeabilities, i.e. from $10^{-8}$ to 1 and permeabilities that contain jumps and anisotropy of similar sizes. Both preconditioners work well; though $\mathcal{B}_{1}$ appears to result in a smaller condition number. $\mathcal{B}_{1}$ has the advantage of performing better, directly, in the case of matrix-valued $K$. However, $\mathcal{B}_{2}$ may be of interest in multi-physics codes where scaling of the $B^{T}B$ component is not desired; in this case, significantly improved condition numbers have been observed when using the exact Schur complement to approximate $-\operatorname{\nabla\cdot}(K\nabla)$ in the presence of matrix-valued $K$. That is, in the presence of anisotropy, a crucial component is a to find a proper local operator representing $B(\alpha A)^{-1}B^{T}$.

7. Acknowledgements

The authors gratefully acknowledge helpful discussions with Ragnar Winther111Dept. of Mathematics, University of Oslo. Oslo, Norway. in addition to productive exchanges with Marie E. Rognes222Dept. of Numerical Anal. and Sci. Comput. , Simula Research Laboratory. Fornebu, Norway and Walter Zulehner333 Johannes Kepler University Linz, Institute of Computational Mathematics Linz, Austria. The work of Miroslav Kuchta was supported by the Research Council of Norway (NFR) grant 280709. The work of Travis Thompson was supported by the Research Council of Norway under the FRINATEK Young Research Talents Programme through project #250731/F20 (Waterscape).