Efficient solvers for hybridized three-field mixed finite element coupled poromechanics

Abstract

We consider a mixed hybrid finite element formulation for coupled poromechanics. A stabilization strategy based on a macro-element approach is advanced to eliminate the spurious pressure modes appearing in undrained/incompressible conditions. The efficient solution of the stabilized mixed hybrid block system is addressed by developing a class of block triangular preconditioners based on a Schur-complement approximation strategy. Robustness, computational efficiency and scalability of the proposed approach are theoretically discussed and tested using challenging benchmark problems on massively parallel architectures.

Introduction

We focus on a three-field mixed (displacement–velocity–pressure) formulation of classical linear poroelasticity [1], [2]. Let $\Omega \subset {\mathbb{R}}^d$ ($d=2,3$) and $\Gamma$ be the domain occupied by the porous medium and its Lipschitz-continuous boundary, respectively, with $\mathit{x}$ the position vector in ${\mathbb{R}}^d$. We denote time with $t$, belonging to an open interval $\mathcal{I}=\left(0,t_{max}\right)$. The boundary is decomposed as $\Gamma =\overline{{\Gamma }_u\cup {\Gamma }_{\sigma }}=\overline{{\Gamma }_p\cup {\Gamma }_q}$, with ${\Gamma }_u\cap {\Gamma }_{\sigma }={\Gamma }_p\cap {\Gamma }_q=0̸$, and $\mathit{n}$ denotes its outer normal vector. Assuming quasi-static, saturated, single-phase flow of a slightly compressible fluid, the set of governing equations consists of a conservation law of linear momentum and a conservation law of mass expressed in mixed form, i.e., introducing Darcy’s velocity as an additional unknown. The strong form of the initial–boundary value problem (IBVP) consists of finding the displacement $\mathit{u}:\overline{\Omega }\times \mathcal{I}\to {\mathbb{R}}^d$, the Darcy velocity $\mathit{q}:\overline{\Omega }\times \mathcal{I}\to {\mathbb{R}}^d$, and the excess pore pressure $p:\overline{\Omega }\times \mathcal{I}\to \mathbb{R}$ that satisfy:

$$\nabla \cdot \left({ℂ}_{dr}:{\nabla }^s\mathit{u}-bp\mathit{1}\right)=\mathit{0}\text{ in }\Omega \times \mathcal{I}\text{(equilibrium)},$$$$\mu {\mathit{\kappa }}^{-1}\cdot \mathit{q}+\nabla p=\mathit{0}\text{ in }\Omega \times \mathcal{I}\text{(Darcy’s law)},$$$$b\nabla \cdot \dot{\mathit{u}}+S_{ϵ}\dot{p}+\nabla \cdot \mathit{q}=f\text{ in }\Omega \times \mathcal{I}\text{(continuity)}.$$

Here, ${ℂ}_{dr}$ is the rank-four elasticity tensor, ${\nabla }^s$ is the symmetric gradient operator, $b$ is the Biot coefficient, and $\mathit{1}$ is the rank-two identity tensor; $\mu$ and $\mathit{\kappa }$ are the fluid viscosity and the rank-two permeability tensor, respectively; $S_{ϵ}$ is the constrained specific storage coefficient, i.e. the reciprocal of Biot’s modulus, and $f$ the fluid source term. The following set of boundary and initial conditions complete the formulation:

$$\mathit{u}=\bar{\mathit{u}}\text{ on }{\Gamma }_u\times \mathcal{I},$$$$\left({ℂ}_{dr}:{\nabla }^s\mathit{u}-bp\mathit{1}\right)\cdot \mathit{n}=\bar{\mathit{t}}\text{ on }{\Gamma }_{\sigma }\times \mathcal{I},$$$$\mathit{q}\cdot \mathit{n}=\bar{q}\text{ on }{\Gamma }_q\times \mathcal{I},$$$$p=\bar{p}\text{ on }{\Gamma }_p\times \mathcal{I},$$$$p(\mathit{x},0)=p_0\mathit{x}\in \overline{\Omega },$$

where $\bar{\mathit{u}}$, $\bar{\mathit{t}}$, $\bar{q}$, and $\bar{p}$ are the prescribed boundary displacements, tractions, Darcy velocity and excess pore pressure, respectively, whereas $p_0$ is the initial excess pore pressure. More precisely, the initial condition should be given as

$$b\nabla \cdot \mathit{u}(\mathit{x},0)+S_{ϵ}p(\mathit{x},0)=b\nabla \cdot {\mathit{u}}_0+S_{ϵ}{p}_0,\mathit{x}\in \overline{\Omega },$$

with ${\mathit{u}}_0$ the initial displacement field, i.e. specifying the initial fluid content of the medium [1]. However, in practical simulations, the pressure is often measured or computed through the hydrostatic assumption, and the initial displacement is then obtained so as to satisfy Eq. (1a)—see, e.g. [3], [4]. We refer the reader to [5] for a rigorous discussion on this issue.

Let us denote with ${\mathit{H}}^1\left(\Omega \right)$ the Sobolev space of vector functions whose first derivatives are square-integrable, i.e., they belong to the Lebesgue space $L^2\left(\Omega \right)$; and let $\mathit{H}(\text{div};\Omega )$ be the Sobolev space of vector functions with square-integrable divergence. Introducing the spaces:

$$\mathcal{U}=\{\mathit{u}\in {\mathit{H}}^1\left(\Omega \right)|\mathit{u}{|}_{{\Gamma }_u}=\bar{\mathit{u}}\},{\mathcal{U}}_0=\{\mathit{u}\in {\mathit{H}}^1\left(\Omega \right)|\mathit{u}{|}_{{\Gamma }_u}=\mathit{0}\},$$$$\mathcal{Q}=\{\mathit{q}\in \mathit{H}(\text{div};\Omega )|\mathit{q}\cdot \mathit{n}{|}_{{\Gamma }_q}=\bar{q}\},{\mathcal{Q}}_0=\{\mathit{q}\in \mathit{H}(\text{div};\Omega )|\mathit{q}\cdot \mathit{n}{|}_{{\Gamma }_q}=0\},$$$$\mathcal{P}=\{p\in L^2\left(\Omega \right)\},$$

the weak form of the IBVP (1) reads: find $\{\mathit{u}\left(t\right),\mathit{q}\left(t\right),p\left(t\right)\}\in \mathcal{U}\times \mathcal{Q}\times \mathcal{P}$ such that $\forall t\in \mathcal{I}$:

$${({\nabla }^s\mathit{\eta },{ℂ}_{\text{dr}}:{\nabla }^s\mathit{u})}_{\Omega }-{(\text{div}\mathit{\eta },bp)}_{\Omega }={(\mathit{\eta },\bar{\mathit{t}})}_{{\Gamma }_{\sigma }}\forall \mathit{\eta }\in {\mathcal{U}}_0,$$$${(\mathit{\varphi },\mu {\mathit{\kappa }}^{-1}\cdot \mathit{q})}_{\Omega }-{(\text{div}\mathit{\varphi },p)}_{\Omega }=-{(\mathit{\varphi }\cdot \mathit{n},\bar{p})}_{{\Gamma }_p}\forall \mathit{\varphi }\in {\mathcal{Q}}_0,$$$${(\chi ,b\text{div}\dot{\mathit{u}})}_{\Omega }+{(\chi ,\text{div}\mathit{q})}_{\Omega }+{(\chi ,S_{ϵ}\dot{p})}_{\Omega }={(\chi ,f)}_{\Omega }\forall \chi \in L^2\left(\Omega \right),$$

where ${(\cdot ,\cdot )}_{\Omega }$ denote the inner products of scalar functions in $L^2\left(\Omega \right)$, vector functions in ${\left[L^2\left(\Omega \right)\right]}^d$, or second-order tensor functions in ${\left[L^2\left(\Omega \right)\right]}^{d\times d}$, as appropriate, and ${(\cdot ,\cdot )}_{{\Gamma }_{\ast }}$ denote the inner products of scalar functions or vector functions on the boundary ${\Gamma }_{\ast }$. For the analysis of the well-posedness of the Biot continuous problem (1) in weak form (5) based on the displacement–velocity–pressure formulation, the reader is referred to [6].

A widely-used discrete version of the weak form (5) is based on low-order elements. Precisely, lowest-order continuous (${\mathbb{Q}}_1$), lowest-order Raviart–Thomas (${\mathbb{RT}}_0$), and piecewise constant (${\mathbb{P}}_0$) spaces are often used for the approximation of displacement, Darcy’s velocity, and fluid pore pressure, respectively. The attractive features of this choice are element-wise mass conservation and robustness with respect to highly heterogeneous hydromechanical properties, such as high-contrast permeability fields typically encountered in real-world applications. Another attractive feature stems from the hybridization of the mixed three-field formulation, as proposed for instance in [7]. The hybridized formulation is obtained by (i) considering one degree of freedom per edge/face per element for the normal component of Darcy’s velocity, and (ii) introducing one Lagrange multiplier on edges/faces in the computational mesh, i.e. an interface pressure, to enforce velocity continuity. The main advantage of the hybrid formulation is that it is amenable to static condensation.

Unfortunately, the ${\mathbb{Q}}_1$-${\mathbb{RT}}_0$-${\mathbb{P}}_0$ discretization spaces do not intrinsically satisfy the LBB condition in the undrained/incompressible limit [6], [8]. This can result in spurious modes for the pressure, with non physical oscillations of the discrete solution. Different stabilization strategies have been proposed in the literature. They can be essentially classified in two groups based on whether they: (1) enrich the discretization spaces to guarantee the LBB condition, or (2) introduce a proper stabilization term to restore the saddle-point problem solvability. In the context of Biot’s poroelasticity, the first strategy is followed for instance in [9], [10], where proper bubble functions are introduced to enlarge the space used for the displacement approximation. This is a mathematically elegant and robust stabilization technique, but it can negatively impact the algebraic structure of the resulting discrete problem, with a possible degradation of the solver computational efficiency. The second approach, used for instance in [11], [12], has a much smaller impact on the algebraic structure of the problem, but depends on the choice of appropriate stabilization coefficients that typically introduce some numerical diffusion. Such coefficients should be properly tuned to guarantee the stabilization effectiveness with no detrimental effect on the solution accuracy. In this work, we adopt the second approach by proposing a local pressure jump stabilization technique based on a macro-element approach [12], [13] that is applicable to both mixed and mixed hybrid formulations.

Then, we concentrate on the efficient solution of the non-symmetric algebraic systems obtained by the application of the stabilized formulation. Several strategies have been already developed for the two-field and mixed three-field formulations, with most of the efforts towards preconditioned Krylov solvers and multigrid methods [6], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26], [27], [28], [29], [30], [31], [32]. Some authors also focused on sequential-implicit approaches [3], [33], [34], [35], [36], [37], [38], [39], [40], [41], [42], where the discrete poromechanical equilibrium equation and the Darcy flow sub-problem are addressed independently, iterating until convergence. Here, a class of block-triangular preconditioners for accelerating the iterative convergence by Krylov subspace methods is proposed for the stabilized mixed hybrid approach. We prove that the hybridization of the classical three-field mixed formulation brings better algebraic properties for the resulting discrete problem, which are exploited by the proposed iterative solver. Performance and robustness of the algorithms are demonstrated in weak and strong scaling studies including both theoretical and field application benchmarks. Finally, a few concluding remarks close the presentation.