Robust preconditioners for a new stabilized discretization of the poroelastic equations

TL;DR

This paper develops and tests robust block preconditioners for a stabilized discretization of poroelastic equations, ensuring efficiency and stability across various parameters in both 2D and 3D problems.

Contribution

It introduces new block preconditioners that are robust with respect to physical and discretization parameters for a stabilized poroelastic discretization.

Findings

Preconditioners are robust across parameters.

Numerical tests confirm effectiveness in 2D and 3D.

Framework improves computational stability.

Abstract

In this paper, we present block preconditioners for a stabilized discretization of the poroelastic equations developed in [45]. The discretization is proved to be well-posed with respect to the physical and discretization parameters, and thus provides a framework to develop preconditioners that are robust with respect to such parameters as well. We construct both norm-equivalent (diagonal) and field-of-value-equivalent (triangular) preconditioners for both the stabilized discretization and a perturbation of the stabilized discretization that leads to a smaller overall problem after static condensation. Numerical tests for both two- and three-dimensional problems confirm the robustness of the block preconditioners with respect to the physical and discretization parameters.

Tables

Table 1. Table 1: Full bubble system. Iteration counts for the block preconditioners on the 2D Mandel problem with varying discretization parameters.

	${ℬ}_D$
	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$	$\frac{1}{64}$	$\frac{1}{128}$
$0.1$	39	40	40	40	38
$0.01$	26	34	39	39	38
$0.001$	23	23	28	34	37
$0.0001$	21	21	21	21	21

Table 2. Table 2: Bubble-eliminated system. Iteration counts for the block preconditioners on the 2D Mandel problem with varying discretization parameters

	${ℬ}_D^E$
	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$	$\frac{1}{64}$	$\frac{1}{128}$
$0.1$	36	40	43	43	42
$0.01$	26	30	37	40	40
$0.001$	32	29	25	31	35
$0.0001$	34	35	31	25	26

Table 3. Table 3: Full bubble system. Iteration counts for the block preconditioners on the 2D Mandel problem with varying physical parameters K 𝐾 K and ν 𝜈 \nu .

	$\nu =0.0$ and varying $K$
	$1$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$
${ℬ}_D$	23	25	35	38	29	19
${ℬ}_L$	7	11	15	17	15	9
${ℬ}_U$	13	16	17	16	15	7
$\widehat{{ℬ}_D}$	35	33	36	38	29	19
$\widehat{{ℬ}_L}$	14	15	16	18	15	10
$\widehat{{ℬ}_U}$	27	22	17	17	15	8

Table 4. Table 4: Bubble-eliminated system. Iteration counts for the block preconditioners on the 2D Mandel problem with varying physical parameters K 𝐾 K and ν 𝜈 \nu .

	$\nu =0.0$ and varying $K$
	$1$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$
${ℬ}_D^E$	36	36	41	42	26	34
${ℬ}_L^E$	17	17	19	21	18	16
${ℬ}_U^E$	23	22	22	21	17	12
${\widehat{{ℬ}_D}}^E$	36	38	41	43	26	34
${\widehat{{ℬ}_L}}^E$	20	20	20	23	18	17
${\widehat{{ℬ}_U}}^E$	27	27	22	21	17	13

Table 5. Table 5: Full bubble system. Iteration counts for the block preconditioners on the 3D footing problem with varying discretization parameters. ( ∗ * means the direct method for solving diagonal blocks is out of memory).

	${ℬ}_D$
	$\frac{1}{4}$	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$
$0.1$	60	65	65	$\ast$
$0.01$	47	57	68	$\ast$
$0.001$	40	42	49	$\ast$
$0.0001$	40	42	42	$\ast$

Table 6. Table 6: Bubble-eliminated system. Iteration counts for the block preconditioners on the 3D footing problem with varying discretization parameters. ( ∗ * means the direct method for solving diagonal blocks is out of memory).

	${ℬ}_D^E$
	$\frac{1}{4}$	$\frac{1}{8}$	$\frac{1}{16}$	$\frac{1}{32}$
$0.1$	61	65	66	$\ast$
$0.01$	54	58	66	$\ast$
$0.001$	58	58	53	$\ast$
$0.0001$	59	61	60	$\ast$

Table 7. Table 7: Full bubble system. Iteration counts for the block preconditioners on the 3D footing problem with varying physical parameters, K 𝐾 K and ν 𝜈 \nu .

	$\nu =0.2$ and varying $K$
	$1$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$
${ℬ}_D$	28	28	49	68	42	35
${ℬ}_L$	20	20	27	37	26	24
${ℬ}_U$	18	18	26	35	21	14
$\widehat{{ℬ}_D}$	28	28	49	68	42	42
$\widehat{{ℬ}_L}$	20	20	28	37	27	25
$\widehat{{ℬ}_U}$	21	21	27	35	22	24

Table 8. Table 8: Bubble-eliminated system. Iteration counts for the block preconditioners on the 3D footing problem with varying physical parameters, K 𝐾 K and ν 𝜈 \nu .

	$\nu =0.2$ and varying $K$
	$1$	${10}^{-2}$	${10}^{-4}$	${10}^{-6}$	${10}^{-8}$	${10}^{-10}$
${ℬ}_D^E$	33	33	51	66	60	61
${ℬ}_L^E$	20	20	29	43	38	35
${ℬ}_U^E$	20	20	29	41	28	18
${\widehat{{ℬ}_D}}^E$	33	33	51	66	60	61
${\widehat{{ℬ}_L}}^E$	22	22	30	43	38	36
${\widehat{{ℬ}_U}}^E$	22	22	29	41	29	29

Table 9. Table 9: Full bubble system. Iteration counts for the block preconditioners on the 3D footing problem with a varying jump in the physical parameter K 𝐾 K .

	$\nu =0.2$ and $k\left(x\right)={10}^{-10}$ for
$k\left(x\right)$ for $x\ge 0.5$	${10}^{-10}$	${10}^{-8}$	${10}^{-6}$	${10}^{-4}$	${10}^{-2}$	$1$
${ℬ}_D$	35	42	84	98	80	80
${ℬ}_L$	24	27	46	56	51	51
${ℬ}_U$	14	20	38	44	39	39
$\widehat{{ℬ}_D}$	42	44	84	98	80	80
$\widehat{{ℬ}_L}$	25	28	46	56	52	51
$\widehat{{ℬ}_U}$	24	22	39	45	44	44

Table 10. Table 10: Bubble-eliminated system. Iteration counts for the block preconditioners on the 3D footing problem with a varying jump in the physical parameter K 𝐾 K .

	$\nu =0.2$ and $k\left(x\right)={10}^{-10}$ for
$k\left(x\right)$ for $x\ge 0.5$	${10}^{-10}$	${10}^{-8}$	${10}^{-6}$	${10}^{-4}$	${10}^{-2}$	$1$
${ℬ}_D^E$	61	62	115	147	131	132
${ℬ}_L^E$	35	39	74	84	77	78
${ℬ}_U^E$	18	27	54	61	56	57
${\widehat{{ℬ}_D}}^E$	61	62	115	147	131	133
${\widehat{{ℬ}_L}}^E$	36	39	74	84	79	79
${\widehat{{ℬ}_U}}^E$	29	29	55	63	61	60

Equations

$$-\mathrm{div}{\mathbf{\sigma }}^{\prime }+\alpha \nabla p=\rho \mathbf{g},\mathrm{in}\Omega ,$$

$${\mathbf{\sigma }}^{\prime }=2\mu \mathbf{\epsilon }(\mathbf{u})+\lambda \mathrm{div}(\mathbf{u})\mathbf{I},\mathrm{in}\Omega ,$$

$$\mathbf{\epsilon }(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^t),\mathrm{in}\Omega ,$$

$$\mathbf{w}=-\frac{1}{{\mu }_f}\mathbf{K}(\nabla p-{\rho }_f\mathbf{g}),\mathrm{in}\Omega ,$$

$$\frac{\partial }{\partial t}\left(\frac{1}{M}p+\alpha \mathrm{div}\mathbf{u}\right)+\mathrm{div}\mathbf{w}=f,\mathrm{in}\Omega ,$$

$$-\mathrm{div}{\mathbf{\sigma }}^{\prime }+\alpha \nabla p=\rho \mathbf{g},\text{mbox}where{\mathbf{\sigma }}^{\prime }=2\mu \mathbf{\epsilon }(\mathbf{u})+\lambda \mathrm{div}(\mathbf{u})\mathbf{I},$$

$$\frac{\partial }{\partial t}\left(\frac{1}{M}p+\alpha \mathrm{div}\mathbf{u}\right)+\mathrm{div}\mathbf{w}=f,$$

$${\mathbf{K}}^{-1}{\mu }_f\mathbf{w}+\nabla p={\rho }_f\mathbf{g}.$$

$$p=0,\text{mbox}forx\in {\overline{\Gamma }}_t,{\mathbf{\sigma }}^{\prime }\mathbf{n}=0,\text{mbox}forx\in {\Gamma }_t,$$

$$\mathbf{u}=0,\text{mbox}forx\in {\overline{\Gamma }}_c,\mathbf{w}\cdot \mathbf{n}=0,\text{mbox}forx\in {\Gamma }_c,$$

$$\left(\frac{1}{M}p+\alpha \mathrm{div}\mathbf{u}\right)(\mathbf{x},0)=0,\mathbf{x}\in \Omega .$$

$$a(\mathbf{u},\mathbf{v})-(\alpha p,\mathrm{div}\mathbf{v})=(\rho \mathbf{g},\mathbf{v}),\forall \mathbf{v}\in \mathbf{V},$$

$$\left(\frac{1}{M}\frac{\partial p}{\partial t},q\right)+\left(\alpha \mathrm{div}\frac{\partial \mathbf{u}}{\partial t},q\right)+(\mathrm{div}\mathbf{w},q)=(f,q),\forall q\in Q,$$

$$({\mathbf{K}}^{-1}{\mu }_f\mathbf{w},\mathbf{r})-(p,\mathrm{div}\mathbf{r})=({\rho }_f\mathbf{g},\mathbf{r}),\forall \mathbf{r}\in \mathbf{W},$$

$$a(\mathbf{u},\mathbf{v})=2\mu {\int }_{\Omega }\mathbf{\epsilon }(\mathbf{u}):\mathbf{\epsilon }(\mathbf{v})+\lambda {\int }_{\Omega }\mathrm{div}\mathbf{u}\mathrm{div}\mathbf{v},$$

$$\mathbf{V}=\{\mathbf{u}\in {\mathbf{H}}^1(\Omega )|\mathbf{u}{|}_{{\overline{\Gamma }}_c}=0\},$$

$$Q=L^2(\Omega ),$$

$$\mathbf{W}=\{\mathbf{w}\in \mathbf{H}(\mathrm{div},\Omega )|(\mathbf{w}\cdot \mathbf{n}){|}_{{\Gamma }_c}=0\},$$

$${\mathbf{V}}_h\subset \mathbf{V},Q_h\subset Q,{\mathbf{W}}_h\subset \mathbf{W}.$$

$$a({\mathbf{u}}_h^m,{\mathbf{v}}_h)-(\alpha p_h^m,\mathrm{div}{\mathbf{v}}_h)=(\rho \mathbf{g},{\mathbf{v}}_h),\forall {\mathbf{v}}_h\in {\mathbf{V}}_h,$$

$$\left(\frac{1}{M}{p}_h^m,q_h\right)+\left(\alpha \mathrm{div}{\mathbf{u}}_h^m,q_h\right)+\tau (\mathrm{div}{\mathbf{w}}_h^m,q_h)=(\tilde{f},q_h),\forall q_h\in Q_h,$$

$$\tau ({\mathbf{K}}^{-1}{\mu }_f{\mathbf{w}}_h^m,{\mathbf{r}}_h)-\tau (p_h^m,\mathrm{div}{\mathbf{r}}_h)=\tau ({\rho }_f\mathbf{g},{\mathbf{r}}_h),\forall {\mathbf{r}}_h\in {\mathbf{W}}_h,$$

$$\mathcal{A}\left(\begin{array}[]{c}{\bm{u}}_{h}^{b}\\ {\bm{u}}_{h}^{l}\\ p_{h}\\ {\bm{w}}_{h}\end{array}\right)={\bm{b}},\ \ \hbox{with}\ \ \mathcal{A}=\left(\begin{array}[]{cccc}A_{bb}&A_{bl}&\alpha B_{b}^{T}&0\\ A_{bl}^{T}&A_{ll}&\alpha B_{l}^{T}&0\\ -\alpha B_{b}&-\alpha B_{l}&\frac{1}{M}M_{p}&-\tau B_{\bm{w}}\\ 0&0&\tau B_{\bm{w}}^{T}&\tau M_{\bm{w}}\end{array}\right),$$

$$a({\mathbf{u}}_h^b,{\mathbf{v}}_h^b)\to A_{bb},a({\mathbf{u}}_h^l,{\mathbf{v}}_h^b)\to A_{bl},a({\mathbf{u}}_h^l,{\mathbf{v}}_h^l)\to A_{ll},$$

$$-(\mathrm{div}{\mathbf{u}}_h^b,q_h)\to B_b,-(\mathrm{div}{\mathbf{u}}_h^l,q_h)\to B_l,-(\mathrm{div}{\mathbf{w}}_h,q_h)\to B_{\mathbf{w}},$$

$$({\mathbf{K}}^{-1}{\mu }_f{\mathbf{w}}_h,{\mathbf{r}}_h)\to M_{\mathbf{w}},\left(p_h,q_h\right)\to M_p,$$

$$a({\mathbf{u}}_h,{\mathbf{v}}_h)\to A_{\mathbf{u}},-(\mathrm{div}{\mathbf{u}}_h,q_h)\to B_{\mathbf{u}},$$

$$\mathcal{A}^{D}=\left(\begin{array}[]{cccc}D_{bb}&A_{bl}&\alpha B_{b}^{T}&0\\ A_{bl}^{T}&A_{ll}&\alpha B_{l}^{T}&0\\ -\alpha B_{b}&-\alpha B_{l}&\frac{1}{M}M_{p}&-\tau B_{\bm{w}}\\ 0&0&\tau B_{\bm{w}}^{T}&\tau M_{\bm{w}}\end{array}\right).$$

$$\mathcal{A}^{E}=\left(\begin{array}[]{ccc}A_{ll}-A_{bl}^{T}D_{bb}^{-1}A_{bl}&\alpha B_{l}^{T}-\alpha A_{bl}^{T}D_{bb}^{-1}B_{b}^{T}&0\\ -\alpha B_{l}+\alpha B_{b}D_{bb}^{-1}A_{bl}&\frac{1}{M}M_{p}+\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{T}&-\tau B_{\bm{w}}\\ 0&\tau B_{\bm{w}}^{T}&\tau M_{\bm{w}}\end{array}\right).$$

$$\left|\left|\left|{\mathbf{x}}_h\right|\right|\right|:={\left[\Vert {\mathbf{u}}_h{\Vert }_{A_{\mathbf{u}}}^2+c_p^{-1}\Vert p_h{\Vert }_{M_p}^2+\tau \Vert {\mathbf{w}}_h{\Vert }_{M_{\mathbf{w}}}^2+{\tau }^2{c}_p\Vert B_{\mathbf{w}}{\mathbf{w}}_h{\Vert }_{M_p^{-1}}^2\right]}^{1/2},$$

$$\underset{0\ne {\mathbf{x}}_h\in {\mathbf{X}}_h}{\sup }\underset{0\ne {\mathbf{y}}_h\in {\mathbf{X}}_h}{\sup }\frac{(\mathcal{A}{\mathbf{x}}_h,{\mathbf{y}}_h)}{\left|\left|\left|{\mathbf{x}}_h\right|\right|\right|\left|\left|\left|{\mathbf{y}}_h\right|\right|\right|}\le \varsigma ,$$

$$\underset{0\ne {\mathbf{y}}_h\in {\mathbf{X}}_h}{\inf }\underset{0\ne {\mathbf{x}}_h\in {\mathbf{X}}_h}{\sup }\frac{(\mathcal{A}{\mathbf{x}}_h,{\mathbf{y}}_h)}{\left|\left|\left|{\mathbf{x}}_h\right|\right|\right|\left|\left|\left|{\mathbf{y}}_h\right|\right|\right|}\ge \gamma ,$$

$$\Vert {\mathbf{u}}^b{\Vert }_{A_{bb}}\le \Vert {\mathbf{u}}^b{\Vert }_{D_{bb}}\le \eta \Vert {\mathbf{u}}^b{\Vert }_{A_{bb}}\text{and}\Vert \mathbf{u}{\Vert }_{A_{\mathbf{u}}}\le \Vert \mathbf{u}{\Vert }_{A_{\mathbf{u}}^D}\le \eta \Vert \mathbf{u}{\Vert }_{A_{\mathbf{u}}},$$

$$\underset{0\ne {\mathbf{x}}_h\in {\mathbf{X}}_h}{\sup }\underset{0\ne {\mathbf{y}}_h\in {\mathbf{X}}_h}{\sup }\frac{({\mathcal{A}}^D{\mathbf{x}}_h,{\mathbf{y}}_h)}{\Vert {\mathbf{x}}_h{\Vert }_{\mathcal{D}}\Vert {\mathbf{y}}_h{\Vert }_{\mathcal{D}}}\le \tilde{\varsigma },$$

$$\underset{0\ne {\mathbf{y}}_h\in {\mathbf{X}}_h}{\inf }\underset{0\ne {\mathbf{x}}_h\in {\mathbf{X}}_h}{\sup }\frac{({\mathcal{A}}^D{\mathbf{x}}_h,{\mathbf{y}}_h)}{\Vert {\mathbf{x}}_h{\Vert }_{\mathcal{D}}\Vert {\mathbf{y}}_h{\Vert }_{\mathcal{D}}}\ge \tilde{\gamma },$$

$$\mathcal{D}=\left(\begin{array}{cccc}{D}_{bb}& A_{bl}& 0& 0\\ A_{bl}^T& A_{ll}& 0& 0\\ 0& 0& \left(\frac{{\alpha }^2}{{\zeta }^2}+\frac{1}{M}\right)M_p& 0\\ 0& 0& 0& \tau M_{\mathbf{w}}+{\tau }^2{c}_p{A}_{\mathbf{w}}\end{array}\right),$$

$$\underset{0\ne {\mathbf{x}}^E\in {\mathbf{X}}_h^E}{\inf }\underset{0\ne {\mathbf{y}}^E\in {\mathbf{X}}_h^E}{\sup }\frac{({\mathcal{A}}^E{\mathbf{x}}^E,{\mathbf{y}}^E)}{\Vert {\mathbf{x}}^E{\Vert }_{{\mathcal{D}}^E}\Vert {\mathbf{y}}^E{\Vert }_{{\mathcal{D}}^E}}\ge {\gamma }^{\ast },$$

$$\underset{0\ne {\mathbf{x}}^E\in {\mathbf{X}}_h^E}{\sup }\underset{0\ne {\mathbf{y}}^E\in {\mathbf{X}}_h^E}{\sup }\frac{({\mathcal{A}}^E{\mathbf{x}}^E,{\mathbf{y}}^E)}{\Vert {\mathbf{x}}^E{\Vert }_{{\mathcal{D}}^E}\Vert {\mathbf{y}}^E{\Vert }_{{\mathcal{D}}^E}}\le \varsigma ,$$

$$\mathcal{D}^{E}=\left(\begin{array}[]{ccc}A_{ll}-A_{bl}^{T}D_{bb}^{-1}A_{bl}&0&0\\ 0&\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{T}+c_{p}^{-1}M_{p}&0\\ 0&0&\tau M_{\bm{w}}+\tau^{2}c_{p}A_{\bm{w}}\end{array}\right),$$

Full text

1 Introduction

In this work, we study the quasi-static Biot model for soil consolidation, where we assume a porous medium to be linearly elastic, homogeneous, isotropic, and saturated by an incompressible Newtonian fluid. According to Biot’s theory [7], the consolidation process satisfies the following system of partial differential equations (PDEs):

[TABLE]

where $\lambda$ and $\mu$ are the Lamé coefficients, $\alpha$ is the Biot-Willis constant, $M$ is the bulk modulus, $\bm{K}$ is the absolute permeability tensor of the porous medium, $\mu_{f}$ is the viscosity of the fluid, ${\bm{I}}$ is the identity tensor, ${\bm{u}}$ is the displacement vector, $p$ is the pore pressure, ${\bm{\sigma}^{\prime}}$ and ${\bm{\epsilon}}$ are the effective stress and strain tensors for the porous medium, and ${\bm{w}}$ is the percolation velocity, or Darcy’s velocity, of the fluid relative to the soil. The right-hand term ${\bm{g}}$ is the density of applied body forces and the source term $f$ represents a forced fluid extraction or injection process. We consider a bounded open subset $\Omega\subset\mathbb{R}^{d}$, $d=2,3$ with regular boundary $\Gamma$.

In many physical applications, the values of some of the parameters described above may vary over orders of magnitude. For instance, in geophysical applications, the permeability can typically range from $10^{-9}$ to $10^{-21}m^{2}$ [35, 49]. Similarly, in biophysical applications such as in the modeling of soft tissue or bone, the permeability can range from $10^{-14}$ to $10^{-16}m^{2}$ [6, 46, 48]. The Poisson ratio, which is the ratio of transverse strain to axial strain, ranges from $0.1$ to $0.5$ in these applications as well. A Poisson ratio of $0.5$ indicates an incompressible material, at which the linear-elastic term becomes positive semi-definite. Due to the variation in relevant values of these physical parameters, it is important to use discretizations that are stable, independently of the parameters. Therefore, in this work we build upon a parameter-robust discretization introduced in [45].

There are several formulations of Biot’s model, and many stable finite-element schemes have been developed for each of them. For instance, in what is called the two-field formulation (displacement and pressure are unknowns), Taylor-Hood elements which satisfy an appropriate inf-sup condition have been used [41, 42, 43]. Unstable finite-element pairs with appropriate stabilization techniques, such as the MINI element have also been developed [44]. Robust block preconditioners for the two-field formulation were studied in [2]. For three-field formulations (displacement, pressure, and Darcy velocity are unknowns), a parameter-independent approach is found in [29]. There, the parameter-robust stability is studied based on a slightly different norm used here, and robust block diagonal preconditioners are proposed. Another stable discretization for the three-field formulation is Crouzeix-Raviart for displacement, lowest order Raviart-Thomas-Nédélec elements for Darcy’s velocity, and piecewise constants for the pressure [31]. Additionally, a different three-field formulation (displacement, fluid pressure, and total pressure are unknowns) was introduced in [35] and a corresponding parameter-robust scheme is studied in the same paper. For a four-field formulation (stress tensor, fluid flux, displacement, and pore-pressure are unknowns), a stable discretization was developed in [34].

In all the cases above, typical discretizations result in a large-scale linear system of equations to solve at each time step. Such linear systems are usually ill-conditioned and difficult to solve in practice. Also, due to their size, iterative solution techniques are usually considered. One approach to solving the coupled poromechanics equations considered here is a sequential method, such as the fixed stress iteration, which consists of first approximating the fluid part and then the geomechanical part. This is then repeated until the solution has converged to within a specified tolerance (see [3, 5, 9, 10, 40] for details). Another approach, is to solve the linear system simultaneously for all unknowns. Examples of this in poromechanics can be found in [2, 12, 22, 24, 25, 37] and the references within. Analysis from [13, 50] indicates that such a fully-implicit method outperforms the convergence rate of the sequential-implicit methods.

Thus, in this work, we take the latter approach and develop robust block preconditioners (e.g. [19, 20, 26]) to accelerate the convergence of Krylov subspace methods solving the full linear system of equations resulting from the discretization of a three-field formulation of Biot’s model. The proposed preconditioners take advantage of the block structure of the discrete model, decoupling the different fields at the preconditioning stage. Such block preconditioning is primarily attractive due to its simplicity, which allows us to focus on the character of the diagonal blocks, and to leverage extensive work on solving simpler problems. For instance, one can take advantage of algebraic multigrid for some of the blocks [11], or auxiliary space decomposition for others [4, 28]. Finally, since we use a stabilized discretization that is well-posed with respect to the physical and discretization parameters [45], we are able to develop robust block preconditioners that efficiently solve the linear systems, independently of such parameters as well.

The rest of the paper is organized as follows. Section 2 reintroduces the three-field formulation and stabilized finite-element discretization considered. Stability of a perturbation to the finite-element discretization is discussed in Section 3. The block preconditioners are then developed in Section 4, presenting both block diagonal and block triangular approaches. Finally, numerical results confirming the robustness and effectiveness of the preconditioners are shown in Section 5, and concluding remarks are made in Section 6.

2 Three-Field Formulation and its Discretization

The focus of this paper is on the three-field formulation, in which Darcy’s velocity, ${\bm{w}}$, is also a primary unknown in addition to the displacement, ${\bm{u}}$, and pressure, $p$. As a result we have the following system of PDEs:

[TABLE]

This system is often subject to the following set of boundary conditions. Though non-homogeneous boundary conditions can also be used, for the sake of simplicity, we consider homogeneous case in this work:

[TABLE]

where ${\bm{n}}$ is the outward unit normal to the boundary, $\overline{\Gamma}=\overline{\Gamma}_{t}\cup\overline{\Gamma}_{c}$; $\Gamma_{t}$ and $\Gamma_{c}$ are open (with respect to $\Gamma$) subsets of $\Gamma$ with nonzero measure. The initial condition at $t=0$ is given by,

[TABLE]

This yields the following mixed formulation for Biot’s three-field consolidation model:

For each $t\in(0,T]$, find $(\bm{u}(t),p(t),\bm{w}(t))\in\bm{V}\times Q\times\bm{W}$ such that

[TABLE]

where,

[TABLE]

corresponds to linear elasticity and $(\cdot,\cdot)$ denotes the standard inner product on $L^{2}(\Omega)$. The function spaces used in the variational form are

[TABLE]

where ${\bm{H}}^{1}(\Omega)$ is the space of square integrable vector-valued functions whose first derivatives are also square integrable, and $\bm{H}(\operatorname{div},\Omega)$ contains the square integrable vector-valued functions with square integrable divergence.

In [45], we developed a stabilized discretization for the three-field formulation described above. Given a partition of the domain, $\Omega$, into $d$-dimensional simplices, $\mathcal{T}_{h}$, we associate a triple of piecewise polynomial, finite-dimensional spaces,

[TABLE]

More specifically, if we choose a piecewise linear continuous finite-element space, $\bm{V}_{h,1}$, enriched with edge/face (2D/3D) bubble functions, $\bm{V}_{b}$, to form $\bm{V}_{h}=\bm{V}_{h,1}\oplus\bm{V}_{b}$ (see [27, pp. 145-149]), a lowest order Raviart-Thomas-Nédélec space (RT0) for $\bm{W}_{h}$, and a piecewise constant space (P0) for $Q_{h}$, Stokes-Biot stable conditions described in Section 3 are satisfied and the formulation is well-posed. Then, using backward Euler as a time discretization on a time interval $(0,t_{\max{}}]$ with constant time-step size $\tau$, the discrete scheme corresponding to the three-field formulation (1)-(3) reads:

Find $(\bm{u}_{h}^{m},p_{h}^{m},\bm{w}_{h}^{m})\in\bm{V}_{h}\times Q_{h}\times\bm{W}_{h}$ such that

[TABLE]

where $(\widetilde{f},q_{h})=\tau(f,q_{h})+\left(\frac{1}{M}p^{m-1}_{h},q_{h}\right)+\left(\alpha\operatorname{div}\bm{u}_{h}^{m-1},q_{h}\right)$, and $(\bm{u}_{h}^{m},p_{h}^{m},\bm{w}_{h}^{m})$ is an approximation to $\left(\bm{u}(\cdot,t_{m}),p(\cdot,t_{m}),\bm{w}(\cdot,t_{m})\right),$ at time $t_{m}=m\tau,\ m=1,2,\ldots.$ The last equation has been scaled by $\tau$ for symmetry. To simplify the notation, we carry out the following stability analysis for a constant time-step size. However, we note that utilizing a variable time-step size leads to analogous results.

Moreover, this discrete variational form can be represented in block matrix form,

[TABLE]

where ${\bm{u}}_{b}$, ${\bm{u}}_{l}$, $p$, and ${\bm{w}}$ are the unknown vectors for the bubble components of the displacement, the piecewise linear components of the displacement, the pressure, and the Darcy velocity, respectively. The blocks in the definition of matrix $\mathcal{A}$ correspond to the following bilinear forms:

[TABLE]

where ${\bm{u}}_{h}={\bm{u}}_{h}^{l}+{\bm{u}}_{h}^{b}$, ${\bm{u}}_{h}^{l}\in\bm{V}_{h,1}$, ${\bm{u}}_{h}^{b}\in\bm{V}_{b}$, and an analogous decomposition for $\bm{v}_{h}$. We further define two matrices for use later,

[TABLE]

such that A_{\bm{u}}=\left(\begin{array}[]{cc}A_{bb}&A_{bl}\\ A_{lb}&A_{ll}\end{array}\right) and $B_{\bm{u}}=\begin{pmatrix}B_{b}&B_{l}\end{pmatrix}$.

A noteworthy result of [45] is that one can replace the enrichment bubble block, $A_{bb}$, in (5) with a spectrally equivalent diagonal matrix, $D_{bb}:=(d+1)\text{diag}(A_{bb})$, resulting in the following linear operator,

[TABLE]

Not only is the resulting operator sparser than the operator in (5), the stabilization term can be eliminated from the operator in a straightforward way (i.e., static condensation), yielding,

[TABLE]

Thus, we obtain an optimal stable discretization with the lowest possible number of degrees of freedom, equivalent to a discretization with P1-RT0-P0 elements, which itself is not stable [45]. While we have reduced the number of degrees of freedom, we note that the sparsity structure of the stiffness matrix has changed as well. The number of non-zeros added to each row depends on the structure of the mesh. The (1,1) block and the (2,1) block increase in non-zeros per row by the number of elements adjacent to a vertex times dimension. In the worst case scenario for the structured grid formed by division of cubes into tetrahedrons, the (1,1) block grows from 37 non-zeros per row to 81 non-zeros per row. The (1,2) block and (2,2) block increase in non-zeros per row by the number of elements adjacent to each element. For the (1,2) block, this doubles the non-zeros per row. The (2,2) block is originally diagonal so the non-zeros per row increases to (spatial) dimension+2. In all cases, the computational cost of multiplication by the modified matrix has the same asymptotic behavior as the mesh size approaches zero.

In [45], it is discussed that due to the spectral equivalence between $A_{bb}$ and $D_{bb}$, the formulation (6) is still well-posed, and remains well-posed independently of the physical and discretization parameters. In the following section (and appendices), we show this in detail, and prove that formulation (7) is also well posed independently of the physical and discretization parameters.

3 Well-Posedness

The well-posedness of the discretized system provides a convenient framework with which to construct block preconditioners. The discrete system using bubble enriched P1-RT0-P0 finite elements (5), which will be referred to as the “full bubble system”, is shown to be well-posed in [45]. However, since (5) is indefinite, the well-posedness of (7) does not simply follow directly. Therefore, in this section, we show well-posedness of (7), as well as (6), which enables block preconditioners for both the full system and the “bubble-eliminated system” to be constructed using the same framework. Since the proofs are quite technical, we include them in the appendices for completeness. First, we give a short overview of the full bubble system case.

To start, for any symmetric positive definite (SPD) matrix $H$, we define the corresponding inner product as $(\bm{u},\bm{v})_{H}:=(H\bm{u},\bm{v})$ and the induced norm as

$\|\bm{v}\|_{H}^{2}:=(\bm{v},\bm{v})_{H}$. In association with the discretized space, $\bm{X}_{h}:=\bm{V}_{h}\times Q_{h}\times\bm{W}_{h}$, we introduce the following weighted norm, for $\bm{x}_{h}=(\bm{u}_{h},p_{h},\bm{w}_{h})^{T}\in\bm{X}_{h}$,

[TABLE]

where $c_{p}:=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}$ with $\zeta:=\sqrt{\lambda+2\mu/d}$, and $d=2$ or $3$ is the dimension of the problem. Under certain conditions (referred to as Stokes-Biot stability [45, Def. 3.1]) on the space $\bm{X}_{h}$, the block matrix form $\mathcal{A}$ defined in (5) is well-posed with respect to the weighted norm (8), i.e., the following continuity and inf-sup condition hold for $\bm{x}_{h}\in\bm{X}_{h}$ and $\bm{y}_{h}=(\bm{v}_{h},q_{h},\bm{r}_{h})^{T}\in\bm{X}_{h}$,

[TABLE]

with constants $\varsigma>0$ and $\gamma>0$ independent of mesh size $h$, time step size $\tau$, and the physical parameters. As mentioned earlier, these conditions are satisfied by our choice of finite-element spaces.

System (6) satisfies similar continuity and inf-sup conditions. From [45], we know that $A_{bb}$ is spectrally equivalent to $D_{bb}$ and $A_{\bm{u}}$ is spectrally equivalent to $A_{\bm{u}}^{D}=\begin{pmatrix}D_{bb}&A_{bl}\\ A_{lb}&A_{ll}\end{pmatrix}$, specifically,

[TABLE]

where the constant $\eta$ depends only on the shape regularity of the mesh. With the above results, we now state the well-posedness of the system given by (6).

Theorem 3.1**.**

If $(\bm{V}_{h},\bm{W}_{h},Q_{h})$ is Stokes-Biot stable, then:

[TABLE]

where,

[TABLE]

*and $A_{\bm{w}}:=B_{\bm{w}}^{T}M_{p}^{-1}B_{\bm{w}}$. *

Remark 3.2**.**

*In general, Theorem 3.1 implies that if we have a well-posed saddle point problem with an SPD first diagonal block, one can replace the first diagonal block by a spectrally equivalent matrix and the resulting saddle point problem is still well-posed. *

The proof of Theorem 3.1 follows from the framework presented in [29]. We have included a proof in Appendix A, as there are details related to the perturbed bilinear form that are not straightforward. Next, we consider the reduced bubble-eliminated formulation, (7). Here, we denote $\bm{X}_{h}^{E}:={\bm{V}}_{1,h}\times Q_{h}\times\bm{W}_{h}$ the discretized finite-element space after bubble elimination.

Theorem 3.3**.**

If the full system (5) is well-posed, satisfying (9) and (10) with respect to the norm (8), then the bubble-eliminated system, (7), satisfies the following inequalities for $\bm{x}^{E}=(\bm{u}_{l},p_{h},\bm{w}_{h})^{T}\in\bm{X}^{E}_{h}$ and $\bm{y}^{E}=(\bm{v}_{l},q_{h},\bm{r}_{h})^{T}\in\bm{X}^{E}_{h}$,

[TABLE]

and,

[TABLE]

where

[TABLE]

with

[TABLE]

*Thus, (7) is well-posed with respect to the weighted norm (16). *

The proof of Theorem 3.3 is technical due to $\mathcal{A}$ being indefinite. Therefore, it is included in Appendix B for the interested reader.

4 Block Preconditioners

Next, we use the properties of the well-posedness to develop block preconditioners for $\mathcal{A}$ and $\mathcal{A}^{E}$. Following the general framework developed in [12, 36, 39, 50, Hong2019], we first consider block diagonal preconditioners (also known as norm-equivalent preconditioners). Then, we discuss block triangular (upper and lower) preconditioners following the framework developed in [32, 36, 38, 47] for Field-of-Value (FOV) equivalent preconditioners. For both cases, we show that the theoretical bounds on their performance remain independent of the discretization and physical parameters of the problem.

4.1 Block Diagonal Preconditioner

Both the full bubble system, (5), and the bubble-eliminated system, (7), are well-posed, satisfying inf-sup conditions, (10) and (14) respectively. Based on the framework proposed in [36, 39], a natural choice for a norm-equivalent preconditioner is the Riesz operator with respect to the inner product corresponding to respective weighted norm (8) or (16).

4.1.1 Full Bubble System

For the full bubble system, the Riesz operator for (8) takes the following block diagonal matrix form:

[TABLE]

In practice, applying the preconditioner $\mathcal{B}_{D}$ involves the action of inverting the diagonal blocks exactly, which is expensive and sometimes infeasible. Therefore, we replace the diagonal blocks by their spectrally-equivalent symmetric and positive definite approximations,

[TABLE]

Here, $S_{\bm{u}}$, $S_{\bm{w}}$, and $S_{p}$ are spectrally equivalent to the action of the inverse of their respective diagonal blocks in $\mathcal{B}_{D}$,

[TABLE]

where the constants $c_{1,\bm{u}}$, $c_{1,p}$, $c_{1,\bm{w}}$, $c_{2,\bm{u}}$, $c_{2,p}$, and $c_{2,\bm{w}}$ are independent of discretization and physical parameters. In practice, $S_{\bm{u}}$ can be defined by standard multigrid methods. For large values of $\tau$ the matrix $\tau M_{\bm{w}}+\tau^{2}\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}A_{\bm{w}}$ is numerically close to singular and requires special preconditioners. With this in mind, $S_{\bm{w}}$ can be defined by either an HX-preconditioner (Auxiliary Space Preconditioner) [28, 33] or multigrid with special block smoothers [4]. In the case of heterogeneous coefficients, specialized approaches such as in [KrausPrecondHdiv2016] can be used. In the full bubble case, $S_{p}$ is obtained by a diagonal scaling ($M_{p}$ is diagonal when using P0 elements). Thus, $S_{p}=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}M_{p}^{-1}$.

4.1.2 Bubble-Eliminated System

In the bubble-eliminated case, the operator for (16) takes the following block diagonal matrix form:

[TABLE]

Here, $A_{\bm{u}}^{E}=A_{ll}-A_{bl}^{T}D_{bb}^{-1}A_{bl}$ and $A_{p}^{E*}=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)M_{p}+\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{T}$.

Again, we replace the diagonal blocks by their spectrally-equivalent symmetric and positive definite approximations,

[TABLE]

Here, $S_{\bm{u}}^{E}$ and $S_{p}^{E}$ are spectrally equivalent to the action of the inverse of their respective diagonal blocks in, ${\mathcal{B}_{D}^{E}}$,

[TABLE]

where the constants $c_{1,\bm{u}}^{E}$, $c_{1,p}^{E}$, $c_{2,\bm{u}}^{E}$, and $c_{2,p}^{E}$ are independent of discretization and physical parameters. In practice, $S_{\bm{u}}^{E}$ and $S_{\bm{w}}$ can be defined similarly as in the full bubble case, and $S_{p}^{E}$ can be defined through standard multigrid methods, as $A_{p}^{E*}$ is equivalent to a Poisson operator.

Since the preconditioners are derived directly from the well-posedness, they too are robust with respect to the physical and discretization parameters of the problem. When applying the preconditioners to the bubble-eliminated formulation, though, a modest degradation in performance compared to the full bubble system is seen. However, robustness with respect to the parameters remains. These properties are demonstrated in the numerical results section.

4.2 Block Triangular Preconditioner

Next, we consider more general preconditioners, in particular, block upper triangular and block lower triangular preconditioners for the linear system, given by $\mathcal{A}$ or $\mathcal{A}^{E}$, following the framework presented in [2, 32, 36, 38, 47] for FOV-equivalent preconditioners. First, we define the notion of FOV equivalence as in [36]. Given a Hilbert space $X$ and its dual $X^{\prime}$, a left preconditioner, $\mathcal{L}:X^{\prime}\rightarrow X$, and a linear operator, $\mathcal{A}:X\rightarrow X^{\prime}$, are FOV-equivalent if, for any $x\in X$,

[TABLE]

In general, $\mathcal{N}:X^{\prime}\rightarrow X$ can be any SPD operator. Here, we choose $\mathcal{N}$ to be a SPD norm-equivalent preconditioner, and $\Sigma$ and $\Upsilon$ are positive constants, with $\Sigma\leq\Upsilon$. Using this definition, we have the following theorem on the convergence rate of preconditioned GMRES for solving $\mathcal{A}\bm{x}=\bm{f}$.

Theorem 4.1**.**

[18, 21]** If $\mathcal{A}$ and $\mathcal{L}$ are FOV-equivalent and $\bm{x}^{m}$ is the $m$-th iteration of the GMRES method preconditioned with $\mathcal{L}$, and $\bm{x}$ is the exact solution, then

[TABLE]

If the constants $\Sigma$ and $\Upsilon$ are independent of physical and discretization parameters, then $\mathcal{L}$ is a uniform left preconditioner for GMRES.

Remark 4.2**.**

Similar arguments apply to right preconditioners for GMRES, which are used in practice. A right preconditioner, $\mathcal{R}:X^{\prime}\rightarrow X$, and linear operator, $\mathcal{A}:X\rightarrow X^{\prime}$ are FOV-equivalent if, for any $x^{\prime}\in X^{\prime}$,

[TABLE]

4.2.1 Full Bubble System

For the three-field formulation, we first consider the block lower triangular preconditioner,

[TABLE]

and the inexact block lower triangular preconditioner,

[TABLE]

Theorem 4.3**.**

Assuming a shape regular mesh and the discretization described above, there exist constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

Theorem 4.4**.**

Assuming the spectral equivalence relations (19) and (21) hold, $\|I-S_{u}A_{\bm{u}}\|_{A_{\bm{u}}}\leq\rho\leq 0.2228$, and $S_{p}=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}M_{p}^{-1}$, then there exists constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

The proofs of the above two theorems turn out to be a special case of the proofs for the bubble-eliminated system (shown below), and thus are omitted here.

Similar arguments can also be applied to block upper triangular preconditioners. We consider the following for $\mathcal{A}$ in (5),

[TABLE]

and the corresponding inexact preconditioner,

[TABLE]

Parameter robustness for the block upper triangular preconditioners is summarized in the following theorems. Again, as these results are special cases of those in the following section, we only state the results here.

Theorem 4.5**.**

Assuming a shape regular mesh and the discretization described above, there exist constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}^{\prime}=\mathcal{B}_{U}^{-1}\bm{x}$ with $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

Theorem 4.6**.**

Assuming the spectral equivalence relations (19) and (21) hold, $\|I-A_{\bm{u}}S_{\bm{u}}\|_{A_{\bm{u}}}\leq\rho\leq 0.2228$, and $S_{p}=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}M_{p}^{-1}$, then there exists constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}^{\prime}=\mathcal{B}_{U}^{-1}\bm{x}$ with $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

4.2.2 Bubble-Eliminated System

For the three-field formulation, we consider the block lower triangular preconditioner,

[TABLE]

and the inexact block lower triangular preconditioner,

[TABLE]

where $B_{\bm{u}}^{E}=B_{l}-B_{b}D_{bb}^{-1}A_{bl}$. For notational convenience, we define

$A_{p}^{E}=\frac{1}{M}M_{p}+\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{T}$ as the pressure block in the bubble-eliminated system.

Lemma 4.7**.**

If the pair of finite-element spaces $\bm{V}_{h}\times Q_{h}$ is Stokes-stable, i.e., they satisfy the following inf-sup condition [27],

[TABLE]

then, in matrix form, we have,

[TABLE]

with $\gamma_{B}=\gamma_{B}^{0}/\sqrt{d}$. Furthermore, from (11), we have,

[TABLE]

Proof 4.8**.**

Here, we use $\bm{v}$ to denote both the finite-element function and its vector representation. Since $\bm{V}_{h}\times Q_{h}$ is Stokes stable, it satisfies the inf-sup condition in (35), where $\gamma^{0}_{B}>0$ is a constant that does not depend on mesh size. Using the fact that $a(\bm{u},\bm{u})\leq(2\mu+d\lambda)(\epsilon(\bm{u}),\epsilon(\bm{u}))$, we have $\|\bm{v}\|_{A_{\bm{u}}}\leq\sqrt{d}\zeta\|\bm{v}\|_{1}$. Then, for any $p\in Q_{h}$,

[TABLE]

Using (38) and

[TABLE]

*(36) is obtained. Equation (37) follows from (11), the spectral equivalence of $A_{\bm{u}}$ and $A_{\bm{u}}^{D}$. *

In order to prove that (33) and (34) satisfy the requirements to be FOV-equivalent preconditioners for the $\mathcal{A}^{E}$ system we need the following relation for the bubble-eliminated system,

[TABLE]

This is established using Lemma 4.7, the first two by two blocks of Equation (51), and a direct computation. With this result, we now show that (33) satisfies the requirements to be an FOV-equivalent preconditioner for $\mathcal{A}^{E}$.

Theorem 4.9**.**

Assuming a shape regular mesh and the discretization described above, there exists constants $\Sigma$ and $\Upsilon$, independent of discretization or physical parameters, such that, for any $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

Proof 4.10**.**

By direct computation and the Cauchy-Schwarz inequality,

[TABLE]

*By the definitions of matrices $A_{p}^{E_{*}}$, $M_{p}$, and $A_{p}^{E}$, we have *

[TABLE]

Then, using (40) on the $\|(B_{\bm{u}}^{E})(A_{\bm{u}}^{E})^{-1}(B_{\bm{u}}^{E})^{T}p\|_{{(A_{p}^{E*})}^{-1}}$ term and (41) on the $\|(A_{p}^{E})p\|_{{(A_{p}^{E*})}^{-1}}$ term, we obtain,

[TABLE]

Observing that, for $d=2,3$, $a(\bm{v},\bm{v})\leq(2\mu+d\lambda)(\epsilon(\bm{v}),\epsilon(\bm{v}))$ for any $\bm{v}$, and a direct computation of the elimination of the bubble, we have,

[TABLE]

Applying (42) to the $\|B_{\bm{u}}^{E}(A_{\bm{u}}^{E})^{-1}(B_{\bm{u}}^{E})^{T}p\|_{M_{p}^{-1}}$ term with $\bm{v}=(A_{\bm{u}}^{E})^{-1}(B_{\bm{u}}^{E})^{T}p$ gives,

[TABLE]

where we use the fact $\frac{\alpha}{\zeta}\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-\frac{1}{2}}<1$. Rewriting the right hand side,

[TABLE]

The above matrix is SPD, meaning that there is a $\sigma>0$ such that

[TABLE]

Using (39), and the definition of $\|p\|_{A_{p}^{E}}$, we get

[TABLE]

*where $\Sigma=\sigma\frac{1}{2}\min\{1,\frac{\gamma_{B}^{2}}{\eta^{2}}\}$. This provides the lower bound for the bubble-eliminated case. The upper bound follows from the continuity of each term and the Cauchy-Schwarz inequality. *

Next, we prove that (34) satisfies the requirements to be an FOV-equivalent preconditioner for the $\mathcal{A}^{E}$ system when the inexact diagonal blocks are solved to sufficient accuracy.

Theorem 4.11**.**

Assuming the spectral equivalence relations (21), (24), and (25) hold, $\|I-S_{u}^{E}A_{\bm{u}}^{E}\|_{A_{\bm{u}}^{E}}\leq\rho$ and $\|I-S_{p}^{E}(A_{p}^{E*})\|_{(A_{p}^{E*})}\leq\beta$, with $\rho>0$ and $\beta>0$ sufficiently small, then there exist constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

Proof 4.12**.**

Assume that $\|I-S_{u}^{E}A_{\bm{u}}^{E}\|_{A_{\bm{u}}^{E}}\leq\rho$ and that $\|I-S_{p}^{E}(A_{p}^{E*})\|_{(A_{p}^{E*})}\leq\beta$. By direct computation,

[TABLE]

*Using $\|I-S_{u}^{E}A_{\bm{u}}^{E}\|_{A_{\bm{u}}^{E}}\leq\rho$ and $\|I-S_{p}^{E}A_{p}^{E*}\|_{A_{p}^{E*}}\leq\beta$, on $\|A_{\bm{u}}^{E}u\|_{S_{u}^{E}}$,

$\|B_{\bm{u}}^{E}(I-S_{u}^{E}A_{\bm{u}}^{E})u\|_{S_{p}^{E}}$, $\|B_{\bm{u}}^{E}S_{u}^{E}(B_{\bm{u}}^{E})^{T}p\|_{S_{p}^{E}}$, and $\|A_{p}^{E}p\|_{S_{p}^{E}}$ allows us to change norms and apply (42), (40), and (41) to these terms in the same way as in the previous proof. Thus,*

[TABLE]

Then, rewriting the right hand side,

[TABLE]

where

[TABLE]

If $\beta$ and $\rho$ are sufficiently small, then the above matrix is SPD, and there is a $\sigma>0$ such that

[TABLE]

*where $\Sigma=\sigma\frac{(1-\rho)(1-\beta)}{2}\min\{1,\frac{\gamma_{B}^{2}}{\eta^{2}}\}$. This provides the lower bound. The upper bound follows from the continuity of each term and the Cauchy-Schwarz inequality. *

Remark 4.13**.**

*Values for $\beta$ and $\rho$ that are sufficiently small can be calculated numerically. For example, if $0<\beta=\rho<0.1291$, then the above matrix is SPD. *

Similar arguments can also be applied to block upper triangular preconditioners. We consider the following upper preconditioner for $\mathcal{A}^{E}$ in (7),

[TABLE]

where, again, $A_{\bm{u}}^{E}=A_{ll}-A_{bl}^{T}D_{bb}^{-1}A_{bl}$, $A_{p}^{E*}=\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)M_{p}+\alpha^{2}B_{b}D_{bb}^{-1}B_{b}^{T}$, and

$B_{\bm{u}}^{E}=B_{l}-B_{b}D_{bb}^{-1}A_{bl}$ when preconditioning the bubble-eliminated case. The corresponding inexact preconditioner is given by:

[TABLE]

Parameter robustness is obtained for the block upper triangular preconditioners using the following theorems. The proofs are similar in concept to the proofs for Theorem 4.9 and 4.11 and are, therefore, omitted.

Theorem 4.14**.**

Assuming a shape regular mesh and the discretization described above, then there exist constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}^{\prime}=(\mathcal{B}_{U}^{E})^{-1}\bm{x}$ with $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

Theorem 4.15**.**

Assuming (21), (24), and (25) hold, $\|I-A_{\bm{u}}^{E}S_{\bm{u}}^{E}\|_{A_{\bm{u}}}\leq\rho$ and $\|I-A_{p}^{E*}S_{p}^{E}\|_{A_{p}^{E*}}\leq\beta$ with $\rho>0$ and $\beta>0$ sufficiently small, there exist constants $\Sigma$ and $\Upsilon$, independent of discretization and physical parameters, such that, for any $\bm{x}^{\prime}=(\mathcal{B}_{U}^{E})^{-1}\bm{x}$ with $\bm{x}=(\bm{u},p,\bm{w})^{T}\neq\bm{0}$,

[TABLE]

This shows that the constructed block preconditioners are robust with respect to the physical and discretization parameters of the bubble-eliminated system, (7).

5 Numerical Results

In this section, we illustrate the convergence benefits obtained using the preconditioners presented above. All test problems were implemented in the HAZmath library [30], which contains routines for finite elements, multilevel solvers, and graph algorithms. The numerical tests were performed on a workstation with an 8-core 3GHz Intel Xeon “Sandy Bridge” CPU and 32 GB of RAM per core.

For each test we use flexible GMRES to solve the linear system obtained from both the bubble-enriched P1-RT0-P0 discretization, $\mathcal{A}$ (5), and the bubble-eliminated discretization, $\mathcal{A}^{E}$ (7). A stopping tolerance of $10^{-8}$ was used for the relative residual of the linear system, measured relative to the norm of the right hand side. For the discretization parameters, tests cover different mesh sizes and different time step sizes. To show robustness with respect to the physical parameters, the permeability, $\bm{K}$, and the Poisson ratio $\nu$ are varied. We also consider a 3D test problem where there are jumps in the permeability. In all test cases we consider a diagonal permeability tensor $\bm{K}=k\bm{I}$.The exact solves for the blocks in $\mathcal{B}_{D}$, $\mathcal{B}_{L}$, and $\mathcal{B}_{U}$ are done using the UMFPACK library [14, 15, 16, 17]. For the inexact blocks, $S_{\bm{u}}$ and $S_{\bm{u}}^{E}$ are inverted using GMRES preconditioned with unsmoothed aggregation AMG in a V-cycle, solved to a relative residual tolerance of $10^{-3}$. The $S_{\bm{w}}$ block is solved using an auxiliary space preconditioned GMRES to a relative residual tolerance of $10^{-3}$ [4, 28, 33]. Using a piecewise constant finite-element space for pressure results in a diagonal matrix for $M_{p}$, so the action of $S_{p}$ is directly computed in the full bubble case. In the bubble-eliminated case, $S_{p}^{E}$ is inverted using GMRES preconditioned with unsmoothed aggregation AMG in a V-cycle, solved to a relative residual tolerance of $10^{-3}$.

5.1 Two-Dimensional Test Problem

First, we consider the Mandel problem in two-dimensions, which models an infinitely long saturated porous slab sandwiched between a top and bottom rigid frictionless plate, and is an important benchmarking tool as the analytical solution is known [1]. At time $t=0$, each plate is loaded with a constant vertical force of magnitude $2F$ per unit length as shown in Figure 1. The analytical solution for pressure is given by

[TABLE]

where $p_{0}=\frac{1}{3a}B(1+\nu_{u})F$, $B=1$ is the Skempton’s coefficient, $\nu_{u}=\frac{3\nu+B(1-2\nu)}{3-B(1-2\nu)}$ is the undrained Poisson ratio, $c$ is the consolidation coefficient given by $c=K(\lambda+2\mu)$, and $\alpha_{n}$ are the positive roots to the nonlinear equation,

[TABLE]

Due to symmetry of the problem we only need to solve in the top right quadrant, defined as $\Omega=(0,1)\times(0,1)$. We cover $\Omega$ with a uniform triangular grid by dividing an $N\times N$ uniform square grid into right triangles, where the mesh spacing is defined by $h=\frac{1}{N}$. For the material properties, $\mu_{f}=1$, $\alpha=1$, and $M=10^{6}$, the Lamé coefficients are computed in terms of the Young modulus, $E=10^{4}$, and the Poisson ratio, $\nu$: $\lambda=\frac{E\nu}{(1-2\nu)(1+\nu)}$ and $\mu=\frac{E}{1+2\nu}$.

Table 1 shows iterations counts for the block preconditioners on the full bubble system for different mesh sizes and time-step sizes. Here, we take one time step using Backward Euler. The physical parameters used in these tests were $\nu=0.0$ and $k=10^{-6}$. We see from the relatively consistent iteration counts that the preconditioned system is robust with respect to the discretization parameters. The block upper and lower triangular preconditioners contain more coupling information than the block diagonal preconditioners, and as a result we see that they preform better than the block diagonal preconditioners.

Similar observations are made for Table 2, which shows iteration counts for the block preconditioners on the bubble-eliminated system for different mesh sizes and time-step sizes. We see that using the bubble-eliminated system results in a slight degradation in performance, but nothing significant. It is also important to note that the performance impact of using the inexact block preconditioners is negligible versus using the exact block preconditioners. This implies that the inexact preconditioners could potentially be solved with less strict tolerance, resulting in more computational efficiency.

Table 3 and Table 4 show iteration counts for the block preconditioners when the physical values of $\nu$ and $K$ are varied for the full bubble system and bubble-eliminated system. The mesh size is fixed to $h=\frac{1}{128}$, and the time-step size is $\tau=0.01$. Again, we observe robustness, this time with respect to the physical parameters. The use of inexact preconditioners and the bubble elimination have minimal impact on the performance. In the limit of impermeability ($k\rightarrow 0$), or in the limit of the Poisson ratio approaching $0.5$, the three-field Biot model limits to the Stokes’ Equation. An interesting result is the better performance when the system is approaching this case.

Finally, Figure 2 shows the time scaling with respect to mesh size for the three different inexact preconditioners for the full-bubble and bubble-eliminated systems. The timings scale on the order of $O(N\log N)$ which is nearly optimal. We also see that while a single iteration of the block lower or block upper triangular preconditioner will take longer than that of a block diagonal iteration, the fewer required iterations of the block triangular preconditioners results in a net savings in total computational time. The bubble-eliminated system, being a smaller system than the full-bubble system, takes less time to solve. Figure 2 shows that solving the bubble-eliminated system is nearly ten times faster than solving the full-bubble system.

5.2 Three-Dimensional Test Problem

Next, we consider a footing problem in three-dimensions as seen in [23]. The domain is a unit cube modeling a block of porous soil. A uniform load, $\bm{\sigma_{0}}$, of intensity $3\times 10^{4}$ per unit area is applied in a square of size $0.5\times 0.5$ in the middle of the top face. The base of the domain is assumed to be fixed while the rest of the domain is free to drain. The material properties used are $\mu_{f}=1$, $\alpha=1$, and $M=10^{6}$, the Lamé coefficients are computed in terms of the Young modulus and the Poisson ratio as in the 2D problem.

Table 5 and Table 6 show iteration counts for the block preconditioners on both systems, while varying the discretization parameters, mesh size and time-step size. Again, one step of Backward Euler is used to test the preconditioners. The physical parameters for these tables were $\nu=0.2$ and $k=10^{-6}$. Here the benefits of the inexact preconditioners becomes clear, as the exact preconditioners could not be used on the two largest meshes due to memory limitations. The iteration counts confirm that the preconditioned system is robust with respect to the discretization parameters even in three dimensions.

Table 7 and Table 8 show the results when the physical parameters are varied. The mesh size is fixed to $h=\frac{1}{16}$, and the time step size is $\tau=0.01$. Again, we see that the preconditioned system is robust with respect to the physical parameters, and that the use of the inexact preconditioners has little impact on the required iterations. The bubble-eliminated system shows performance that is overall similar to the full bubble system.

Similarly to Figure 2, Figure 4 shows time scaling with respect to mesh size for the three different inexact preconditioners for the full-bubble and bubble-eliminated systems, again showing a nearly optimal scaling of $O(N\log N)$. The time comparison between the three different inexact preconditioners again demonstrate that the block lower and block upper triangular preconditioners are faster than the block diagonal preconditioner despite being more expensive per iteration. Finally, we see that solving the bubble-eliminated system is faster than solving the full-bubble system as expected.

In order to show the full capabilities of the preconditioners, we test on the 3D footing problem when there is a spatially-dependent jump in the value for the permeability tensor $\bm{K}$. The permeability tensor, $\bm{K}=k(x)\bm{I}$, is defined so that $k(x)=10^{-10}$ when $x<0.5$ and varied for $x\geq 0.5$. The results are shown in Table 9 and Table 10. The Poisson ratio is $\nu=0.2$, the mesh size is fixed to $h=\frac{1}{16}$, and the time-step size is $\tau=0.01$. Note that the size of the jump increases from left to right in the table. We see that, at the beginning, the iteration counts for the preconditioned system increases when the jump gets larger. However, it stabilizes as the jump gets larger and, more importantly, the iterations are bounded from above. This is consistent with our theoretical results that there is an upper bound on the condition number or field-of-values for the preconditioned system.

6 Conclusions

The stability and well-posedness of the discrete problem provides a foundation for designing robust preconditioners. Thus, we are able to develop block preconditioners which yield uniform convergence rates for GMRES. These preconditioners are robust with respect to both the physical and the discretization parameters, making it attractive for problems in poromechanics, such as Biot’s consolidation model considered here. Moreover, the bubble-eliminated system has the same number of degrees of freedom as a P1-RT0-P0 discretized system, yet it is well-posed independent of the physical and discretization parameters, and it attains performance similar to the full bubble enriched P1-RT0-P0 system. Due to the lower number of degrees of freedom, though, the solution time is faster than the fully-stabilized system.

Future work involves developing monolithic multigrid methods for the stabilized discretization of the three-field Biot model presented in [45]. The block preconditioners presented here can then be used as a relaxation step in the monolithic multigrid method, and the overall performance will be compared against this work as stand alone preconditioners. Additionally, other test problems including systems with fractures or other nonlinear behavior will be considered.

Appendix A Proof of Theorem 3.1

The following lemmas are useful for the following proofs.

Lemma A.1**.**

Given the system defined in (5),

[TABLE]

*where $(-\operatorname{div}\bm{v},q)\rightarrow B$, with $q\in Q_{h}$ and $\bm{v}\in\bm{V}_{h}$. *

Proof A.2**.**

By direct computation,

[TABLE]

*where $P_{Q_{h}}$ is the $L^{2}$-projection from $Q$ onto $Q_{h}$. As an abuse of notation, we use $\bm{v}$ for the corresponding vector representation and write the above inequality in matrix form, concluding that *

[TABLE]

Corollary A.3**.**

Considering only the bubble component for Lemma A.1, we have

[TABLE]

Proof A.4**.**

*The first inequality follows the same arguments as the proof of Lemma A.1 and the second inequality comes from the spectral equivalence of $A_{bb}$ and $D_{bb}$, i.e., (11). *

With the above result, we now show the well-posedness of the system given by Theorem 3.1, restated here.

Theorem A.5**.**

If $(\bm{V}_{h},\bm{W}_{h},Q_{h})$ is Stokes-Biot stable, then:

[TABLE]

where,

[TABLE]

*and $A_{\bm{w}}:=B_{\bm{w}}^{T}M_{p}^{-1}B_{\bm{w}}$. *

Proof A.6**.**

*From (38), for a given $p\in Q_{h}$, there exists $\bm{z}\in\bm{V}_{h}$, such that

$(p,B_{\bm{u}}\bm{z})\geq\frac{\gamma_{B}}{\eta\zeta}\|p\|_{M_{p}}^{2}$ and $\|\bm{z}\|_{A_{\bm{u}}^{D}}=\|p\|_{M_{p}}$. Let $\bm{v}=\bm{u}-\theta_{1}\bm{z}$, $\bm{r}=\bm{w}$, and

$q=-p-\theta_{2}\tau\operatorname{div}\bm{w}$, then, by the Cauchy-Shwarz and Young’s inequality,*

[TABLE]

Combining terms and applying Lemma A.1,

[TABLE]

Choosing $\theta_{1}=\frac{\alpha\gamma_{B}}{2\eta\zeta}$ and $\theta_{2}=\frac{1}{2}\left(\frac{\alpha^{2}}{\zeta^{2}}+\frac{1}{M}\right)^{-1}$,

[TABLE]

*where $\tilde{\gamma}=\min\left\{\frac{1}{6},\frac{3\gamma_{B}^{2}}{8\eta^{2}}\right\}$. *

Appendix B Proof of Theorem 3.3

To start, we use a result from [8], which is restated here for convenience.

Proposition B.1** (Proposition 3.4.5 in [8]).**

Let $B$ be an $m\times n$ matrix, $S_{X}$ be an $n\times n$ symmetric positive definite matrix, and $S_{Y}$ be an $m\times m$ symmetric positive definite matrix. Define the following norms: $\|\bm{x}\|_{X}^{2}:=(S_{X}\bm{x})^{T}(S_{X}\bm{x})$ and $\|\bm{y}\|_{Y}^{2}:=(S_{Y}\bm{y})^{T}(S_{Y}\bm{y})$ for $\bm{x}\in\mathbb{R}^{n},\ \bm{y}\in\mathbb{R}^{m}$, and let $\beta$ be defined as

[TABLE]

*where $K:=\ker B$ and $H:=\ker B^{T}$. Then, $\beta$ coincides with the smallest positive singular value of the matrix $S_{Y}BS_{X}$. *

With the above result, we now show the well-posedness of the bubble-eliminated system given by Theorem 3.3, restated here.

Theorem B.2**.**

If the full system (5) is well-posed, satisfying (9) and (10) with respect to the norm (8), then the bubble-eliminated system, (7), satisfies the following inequalities for $\bm{x}^{E}=(\bm{u}_{l},p_{h},\bm{w}_{h})^{T}\in\bm{X}^{E}_{h}$ and $\bm{y}^{E}=(\bm{v}_{l},q_{h},\bm{r}_{h})^{T}\in\bm{X}^{E}_{h}$,

[TABLE]

and,

[TABLE]

where,

[TABLE]

with

[TABLE]

*Thus, (7) is well-posed with respect to the weighted norm (16). *

Proof B.3**.**

The matrix $\mathcal{A}^{D}$ given in (6) affords the following decomposition,

[TABLE]

with

[TABLE]

and

[TABLE]

Note that $\mathcal{A}^{E}$ is a sub-matrix of $\mathcal{S}$. Then, looking at the inf-sup condition (13) . we have, for $\bm{x}=(\bm{u}_{b},\bm{u}_{l},p_{h},\bm{w}_{h})^{T}\in\bm{X}_{h}$ and $\bm{y}=(\bm{v}_{b},\bm{v}_{l},q_{h},\bm{r}_{h})^{T}\in\bm{X}_{h}$,

[TABLE]

where $\bm{\xi}=\tilde{L}^{T}\bm{x}$, $\bm{\varphi}=L^{T}\bm{y}$. We will proceed by showing that $L^{-1}\mathcal{D}L^{-T}$ and $\tilde{L}^{-1}\mathcal{D}\tilde{L}^{-T}$ are spectrally equivalent to the following block diagonal matrix,

[TABLE]

Note that $\mathcal{D}^{E}$, corresponding to the weighted norm on the bubble-eliminated system (16), is a sub-matrix of $\tilde{\mathcal{D}}$.

By direct computation, the Cauchy-Schwarz inequality, Young’s inequality, and use of Corollary A.3 we have, for $\bm{x}=(\bm{u}_{b},\bm{u}_{l},p_{h},\bm{w}_{h})^{T}\in\bm{X}_{h}$,

[TABLE]

*where $A_{\bm{u}}^{E}=A_{ll}-A_{bl}^{T}D_{bb}^{-1}A_{bl}$. Thus, we get that *

[TABLE]

Similarly,

[TABLE]

yielding,

[TABLE]

With (52) and (53) we have that $L^{-1}\mathcal{D}L^{-T}$ is spectrally equivalent to $\tilde{\mathcal{D}}$. For the $\tilde{L}^{-1}\mathcal{D}\tilde{L}^{-T}$ operator, by direct computation and the Cauchy-Schwarz inequality, we have

[TABLE]

and the rest of the proof for the lower bound follows exactly as it does in the $L^{-1}\mathcal{D}L^{-T}$ case. Similarly for the upper bound, we have

[TABLE]

and the rest of the proof for the upper bound follows from the $L^{-1}\mathcal{D}L^{-T}$ case. Thus $L^{-1}\mathcal{D}L^{-T}$ and $\tilde{L}^{-1}\mathcal{D}\tilde{L}^{-T}$ are spectrally equivalent to the block diagonal matrix $\tilde{\mathcal{D}}$. We then write, $\forall\bm{x},\bm{y}$,

[TABLE]

Since the maps $\tilde{L}^{T}:\bm{x}\mapsto\bm{\xi}$ and $L^{T}:\bm{y}\mapsto\bm{\varphi}$ are one-to-one,

[TABLE]

where $\gamma^{*}=\frac{\tilde{\gamma}}{16}$.

Evoking Proposition B.1 (Proposition 3.4.5 in [8]), we know that the smallest singular value of $\tilde{\mathcal{D}}^{-1/2}\mathcal{S}\tilde{\mathcal{D}}^{-1/2}$ is bounded from below by a fixed positive constant. The matrix,

[TABLE]

is a block diagonal matrix with $(\mathcal{D}^{E})^{-1/2}\mathcal{A}^{E}(\mathcal{D}^{E})^{-1/2}$ as a submatrix on the diagonal. Then, the smallest singular value of $(\mathcal{D}^{E})^{-1/2}\mathcal{A}^{E}(\mathcal{D}^{E})^{-1/2}$ is bounded from below by a fixed positive constant. Therefore, we arrive at equation (14), for $\bm{x}^{E}=(\bm{u}_{l},p_{h},\bm{w}_{h})^{T}\in\bm{X}^{E}_{h}$ and $\bm{y}^{E}=(\bm{v}_{l},q_{h},\bm{r}_{h})^{T}\in\bm{X}^{E}_{h}$,

[TABLE]

The upper bound follows from the following set of inequalities,

[TABLE]

*which results in (15). *