Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems

TL;DR

This paper introduces adaptive asynchronous time-stepping and stopping criteria for fixed-stress iterative schemes in coupled poromechanics, utilizing a posteriori error estimates to improve accuracy and efficiency in simulations.

Contribution

It develops a posteriori error estimates and adaptive algorithms for fixed-stress schemes, enabling more efficient and accurate coupled poromechanics simulations.

Findings

Guaranteed and fully computable error bounds are derived.

Adaptive algorithms improve convergence and computational efficiency.

Numerical experiments confirm the effectiveness of the proposed methods.

Abstract

In this paper we develop adaptive iterative coupling schemes for the Biot system modeling coupled poromechanics problems. We particularly consider the space-time formulation of the fixed-stress iterative scheme, in which we first solve the problem of flow over the whole space-time interval, then exploiting the space-time information for solving the mechanics. Two common discretizations of this algorithm are then introduced based on two coupled mixed finite element methods in-space and the backward Euler scheme in-time. Therefrom, adaptive fixed-stress algorithms are build on conforming reconstructions of the pressure and displacement together with equilibrated flux and stresses reconstructions. These ingredients are used to derive a posteriori error estimates for the fixed-stress algorithms, distinguishing the different error components, namely the spatial discretization, the temporal discretization, and the fixed-stress iteration components. Precisely, at the iteration $k\geq 1$ of the adaptive algorithm, we prove that our estimate gives a guaranteed and fully computable upper bound on the energy-type error measuring the difference between the exact and approximate pressure and displacement. These error components are efficiently used to design adaptive asynchronous time-stepping and adaptive stopping criteria for the fixed-stress algorithms. Numerical experiments illustrate the efficiency of our estimates and the performance of the adaptive iterative coupling algorithms.

Figures

Tables

Table 1. Table 1 : The three computations in test problem 2.

Algorithm	adaptive		standard
User-weights	${\gamma }_{\text{tm,J}}=0.9$, ${\Gamma }_{\text{tm,J}}=1.1$	${\gamma }_{\text{tm,J}}=0.5$, ${\Gamma }_{\text{tm,J}}=1.5$	none
Tolerance	${\gamma }_{\text{it}}=0.5$	${\gamma }_{\text{it}}=0.5$	${\text{err}}_{\text{FS}}^k\le {10}^{-5}$
Nb. iterations	$\approx 30$	$\approx 30$	$\approx 136$
Nb. unknowns	571201	864341	1920375
Tot. estimate	0.7943	0.723	0.638

Equations

$$-\nabla \cdot (\mathbf{\theta }(\mathbf{u})-\alpha p\mathbf{I})$$

$${\partial }_t\phi (p,\mathbf{u})-\nabla \cdot (\mathbf{K}\nabla p)$$

$$\mathbf{u}(\cdot ,0)={\mathbf{u}}_0,p(\cdot ,0)$$

$$\mathbf{u}=0,p$$

$$||\mathbf{v}|{|}_{\text{div},\Omega }:={\left(||\mathbf{v}|{|}^2+||\nabla \cdot \mathbf{v}|{|}^2\right)}^{\frac{1}{2}}.$$

$$Q$$

$$\mathcal{A}\mathbf{\tau }:=\frac{1}{2\mu }\left(\mathbf{\tau }-\frac{\lambda }{2\mu +d\lambda }\text{tr}(\mathbf{\tau })\mathbf{I}\right),$$

$${\int }_0^T\{(c_0+c_r)({\partial }_{t}p,q)+\frac{c_r}{d\alpha }({\partial }_t\mathbf{\sigma },q\mathbf{I})+(\nabla \cdot \mathbf{w},q)\}\text{d}t$$

$${\int }_0^T\{({\mathbf{K}}^{-1}\mathbf{w},\mathbf{v})-(p,\nabla \cdot \mathbf{v})\}\text{d}t$$

$$-{\int }_0^T\{(\mathcal{A}\mathbf{\sigma },\mathbf{\tau })+(\mathbf{u},\nabla \cdot \mathbf{\tau })+(\mathbf{\zeta },\mathbf{\tau })+\frac{c_r}{d\alpha }(p\mathbf{I},\mathbf{\tau })\}\text{d}t$$

$${\int }_0^T\{(\nabla \cdot \mathbf{\sigma },\mathbf{z})+(\mathbf{\sigma },\mathbf{\gamma })\}\text{d}t$$

$${\int }_0^T\{(c_0+c_r)({\partial }_t{p}^k,q)+(\nabla \cdot {\mathbf{w}}^k,q)\}\text{d}t$$

$${\int }_0^T\{({\mathbf{K}}^{-1}{\mathbf{w}}^k,\mathbf{v})-(p^k,\nabla \cdot \mathbf{v})\}\text{d}t$$

$$-{\int }_0^T\{(\mathcal{A}{\mathbf{\sigma }}^k,\mathbf{\tau })+({\mathbf{u}}^k,\nabla \cdot \mathbf{\tau })+({\mathbf{\zeta }}^k,\mathbf{\tau })\}\text{d}t$$

$${\int }_0^T\{(\nabla \cdot {\mathbf{\sigma }}^k,\mathbf{z})+({\mathbf{\sigma }}^k,\mathbf{\gamma })\}\text{d}t$$

$$Q_h:=\{q_h\in L^2(\Omega );\forall K\in {\mathcal{T}}_h,q_h{|}_K\in \text{amsmathbb}{P}_0(K)\},$$

$${\mathbf{W}}_h:=\{{\mathbf{v}}_h\in \mathbf{H}(\text{div},\Omega );\forall K\in {\mathcal{T}}_h,{\mathbf{v}}_h{|}_K\in {\mathbf{RTN}}_0(K)\},$$

$${\mathbf{Q}}_h:=\{{\mathbf{z}}_h\in {\mathbf{L}}^2(\Omega );\forall K\in {\mathcal{T}}_h,{\mathbf{z}}_h{|}_K\in {\left[\text{amsmathbb}{P}_0(K)\right]}^d\},$$

$$\text{amsmathbb}{W}_h:=\{{\mathbf{\tau }}_h\in \text{amsmathbb}H(\text{div},\Omega );\forall K\in {\mathcal{T}}_h,{\mathbf{\tau }}_h{|}_K\in {\left[\text{amsmathbb}{P}_1(K)\right]}^{d\times d}\},$$

$$\text{amsmathbb}{Q}_{\text{sk},h}:=\{{\gamma }_h\in [L^2(\Omega ){]}_{\text{sk}}^{d\times d};\forall K\in {\mathcal{T}}_h,{\gamma }_h{|}_K\in [\text{amsmathbb}{P}_0(K){]}_{\text{sk}}^{d\times d}\},$$

$${\partial }_t{v}_h^n:={\partial }_t{v}_{h\tau }{|}_{I_{\text{f}}^n}=\frac{v_h^n-v_h^{n-1}}{{\tau }_{\text{f}}^n}.$$

$$({\mathbf{K}}^{-1}{\mathbf{w}}_h^{k,n},\mathbf{v})-(p_h^{k,n},\nabla \cdot \mathbf{v})=0,$$

$$(c_0+c_r+\beta )({\partial }_t{p}_h^{k,n},q)+(\nabla \cdot {\mathbf{w}}_h^{k,n},q)=(g^n,q)$$

$$+\beta ({\partial }_t{p}_h^{k-1,n},q)-\frac{c_r}{d\alpha }({\partial }_t{\mathbf{\sigma }}_h^{k-1,n},q\mathbf{I}),$$

$$(\mathcal{A}{\mathbf{\sigma }}_h^{k,n},\mathbf{\tau })+({\mathbf{u}}_h^{k,n},\nabla \cdot \mathbf{\tau })+({\mathbf{\zeta }}_h^{k,n},\mathbf{\tau })=-\frac{c_r}{d\alpha }(p_h^{k,n}\mathbf{I},\mathbf{\tau }),$$

$$(\nabla \cdot {\mathbf{\sigma }}_h^{k,n},\mathbf{z})+({\mathbf{\sigma }}_h^{k,n},\mathbf{\gamma })=-({\mathbf{f}}^n,\mathbf{z}),$$

$$({\mathbf{K}}^{-1}{\mathbf{w}}_h^{k,\ell +m},\mathbf{v})-(p_h^{k,\ell +m},\nabla \cdot \mathbf{v})=0,$$

$$(c_0+c_r+\beta )({\partial }_t{p}_h^{k,\ell +m},q)+(\nabla \cdot {\mathbf{w}}_h^{k,\ell +m},q)=(g^{\ell +m},q)$$

$$+\beta ({\partial }_t{p}_h^{k-1,\ell +m},q)-\frac{c_r}{d\alpha }(\frac{{\mathbf{\sigma }}_h^{k-1,\ell +{\delta }_{\text{fm}}}-{\mathbf{\sigma }}_h^{k-1,\ell }}{{\tau }_{\text{m}}^{\ell +{\delta }_{\text{fm}}}},q\mathbf{I}),$$

$$(\mathcal{A}{\mathbf{\sigma }}_h^{k,\ell +{\delta }_{\text{fm}}},\mathbf{\tau })+({\mathbf{u}}_h^{k,\ell +{\delta }_{\text{fm}}},\nabla \cdot \mathbf{\tau })+({\mathbf{\zeta }}_h^{k,\ell +{\delta }_{\text{fm}}},\mathbf{\tau })=-\frac{c_r}{d\alpha }(p_h^{k,\ell +{\delta }_{\text{fm}}}\mathbf{I},\mathbf{\tau }),\forall \mathbf{\tau }\in \text{amsmathbb}{W}_h,$$

$$(\nabla \cdot {\mathbf{\sigma }}_h^{k,\ell +{\delta }_{\text{fm}}},\mathbf{z})+({\mathbf{\sigma }}_h^{k,\ell +{\delta }_{\text{fm}}},\mathbf{\gamma })=-({\mathbf{f}}^{\ell +{\delta }_{\text{fm}}},\mathbf{z}),\forall (\mathbf{z},\mathbf{\gamma })\in {\mathbf{Q}}_h\times \text{amsmathbb}{Q}_{\text{sk},h},$$

$${\mathbf{u}}_h^{k,\ell +m}$$

$${\mathbf{\sigma }}_h^{k,\ell +m}$$

$$-\mathbf{K}\nabla {\tilde{p}}_h^{k,n}$$

$$({\tilde{p}}_h^{k,n},1{)}_K$$

$$\nabla {\tilde{\mathbf{u}}}_h^{k,n}-\frac{c_r}{d\alpha }{p}_h^{k,n}\mathbf{I}-{\mathbf{\zeta }}_h^{k,n}$$

$$\frac{({\tilde{\mathbf{u}}}_h^{k,n},{\mathbf{e}}_i{)}_K}{|K|}$$

$${\tilde{p}}_{h\tau }^k(\cdot ,0)=p_0,\text{ and }{\tilde{\mathbf{u}}}_{h\tau }^k(\cdot ,0)={\mathbf{u}}_0.$$

$${\tilde{p}}_{h\tau }^k(\cdot ,t_{\text{f}}^n)={\tilde{p}}_h^{k,n},{\tilde{\mathbf{u}}}_{h\tau }^k(\cdot ,t_{\text{f}}^n)={\tilde{\mathbf{u}}}_h^{k,n},0\le n\le N_{\text{f}}.$$

$${\gamma }_{\text{tm,J}}{\eta }_{\text{sp,J}}^{k,n}\le {\eta }_{\text{tm,J}}^{k,n}\le {\Gamma }_{\text{tm,J}}{\eta }_{\text{sp,J}}^{k,n},\text{J=P,U},$$

Full text

Key words: Biot’s poro-elasticity problem; mixed finite element method; fixed-stress iterative coupling; space-time scheme; multi-rate scheme; Arnold–Falk–Winther elements; a posteriori error analysis; energy-type estimates; adaptive asynchronous time-stepping; adaptive stopping criteria.

1 Introduction

Let $\Omega$ be an open, bounded and connected domain in ${\amsmathbb R}^{d}$, $d=2,3$, which is assumed to be polygonal with Lipschitz-continuous boundary $\partial\Omega$, and let $T$ be the final simulation time. We consider in this paper the problem of flow in deformable porous media modelled by the quasi-static Biot system [42]: find a displacement ${\mathbf{u}}$ and a pressure $p$ satisfying:

[TABLE]

where $\mathbf{f}$ is the body force and $g$ is the volumetric source term. The function ${\boldsymbol{\theta}}$ denotes the effective stress tensor, i.e., ${\boldsymbol{\theta}}({\mathbf{u}}):=2\mu\boldsymbol{\epsilon}({\mathbf{u}})+\lambda\textnormal{tr}(\boldsymbol{\epsilon}({\mathbf{u}})){\mathbf{I}}$, with $\boldsymbol{\epsilon}({\mathbf{u}})$ is the linearized strain tensor given by $\boldsymbol{\epsilon}({\mathbf{u}}):=(\nabla{\mathbf{u}}+\nabla^{\textnormal{T}}{\mathbf{u}})/2$ and the operator “tr”  denotes the trace of matrices. The coefficients $\mu$ and $\lambda$ are the Lamé parameters, supposed strictly positive constants. The function $\varphi$ denotes the fluid content, i.e., $\varphi(p,{\mathbf{u}}):=c_{0}p+\alpha\nabla{\cdot}{\mathbf{u}}$, where $c_{0}>0$ is the constrained-specific storage coefficient, and $\alpha>0$ is the Biot–Willis constant. The parameter $\mathbf{K}$ is the permeability tensor divided by fluid viscosity; it is a symmetric, bounded, and uniformly positive definite tensor whose terms are for simplicity supposed piecewise constant on the mesh $\mathcal{T}_{h}$ of $\Omega$ defined below and constant in time. Finally, $p_{0}$ is the initial pressure and ${\mathbf{u}}_{0}$ is the initial displacement.

Throughout, we will use the convention that if $V$ is a space of functions, then we designate by $\mathbf{V}$ a space of vector functions having each component in $V$, and we designate by $\amsmathbb{V}$ the space of tensor functions having each component in $V$. Let $D\subset{\amsmathbb R}^{d}$; the space $L^{2}(D)$ is endowed with its natural inner product written $\left(\cdot,\cdot\right)_{D}$ with associated norm denoted by $||\cdot||_{D}$. When the domain $D$ coincides with $\Omega$, the subscript $\Omega$ is dropped. Let $|D|$ be the Lebesgue measure of $D$. We designate by $H^{1}(\Omega)$ the usual Sobolev space and by $H^{1}_{0}(\Omega)$ for its zero-trace subspace. Its norm and semi-norm are written $||\cdot||_{H_{0}^{1}(\Omega)}$ and $|\cdot|_{H_{0}^{1}(\Omega)}$ respectively. In particular, $H^{-1}(\Omega)$ is the dual of $H^{1}_{0}(\Omega)$. Further, let $\mathbf{H}({\textnormal{div}},\Omega)$ be the space of vector-valued functions from $\mathbf{L}^{2}(\Omega)$ that admit a weak divergence in $L^{2}(\Omega)$. Its natural norm is

[TABLE]

We also define $\amsmathbb{H}({\textnormal{div}},\Omega)$ to be the space of tensor-valued functions from $\amsmathbb{L}^{2}(\Omega)$ that admit a weak divergence (by rows) in $\mathbf{L}^{2}(\Omega)$. Then, we set

[TABLE]

To give the mixed formulation of (1.1), we introduce the total stress tensor; ${\boldsymbol{\sigma}}(p,{\mathbf{u}}):={\boldsymbol{\theta}}({\mathbf{u}})-\alpha p{\mathbf{I}}$, and the Darcy velocity; ${\mathbf{w}}:=-\mathbf{K}\nabla p$. Let $c_{r}:=\dfrac{d\alpha^{2}}{2\mu+d\lambda}$, then introduce the fourth-order compliance tensor $\mathcal{A}$ given by

[TABLE]

which is known to be bounded and symmetric definite uniformly with respect to ${\mathbf{x}}\in\Omega$, we can rewrite equations (1.1) in mixed weak sense:

Definition 1.1** (The five-field formulation [2]).**

Assume $\mathbf{f}\in L^{2}(0,T;\left[L^{2}(\Omega)\right]^{d})$, $g\in L^{2}(0,T;L^{2}(\Omega))$, $p_{0}\in H^{1}_{0}(\Omega)$ and ${\mathbf{u}}_{0}\in\mathbf{H}^{1}_{0}(\Omega)$. The fully mixed formulation of (1.1) reads: find $({\boldsymbol{\sigma}},{\mathbf{u}},{\mathbf{w}},p,{\boldsymbol{\zeta}})\in H^{1}(0,T;\amsmathbb{W})\times L^{2}(0,T;\mathbf{Q})\times L^{2}(0,T;\mathbf{W})\times H^{1}(0,T;Q)\times L^{2}(0,T;\amsmathbb{Q}_{\textnormal{sk}})$ such that

[TABLE]

together with the initial condition (1.1d).

The well-posedness and regularity analysis of the Biot equations (1.1) have been addressed in [42]. That of the existence and uniqueness of a weak solution of problem (1.3) have been addressed in [3] (see [2] for more details). Therein, two mixed formulations are discretized with the backward Euler scheme in-time and in-space with mixed finite elements methods, then a posteriori error estimates for their solutions are derived. The main issue arising when apply MFE methods for this problem is that it results a very large system to be solved at each time step (see [3, 15, 45]). This issue together with the fact that flow and mechanics effects act at different time scales, encourage the development of efficient techniques for the resolution of these coupled systems. Splitting-based iterative methods [17, 27, 31, 32] provide one such approach. They adopt the “divide and conquer”  strategy and split the two systems. Then, a sequential approach is used, in that either the problem of flow or the mechanics is solved first followed by solving the other system using the already calculated information, leading to recover the original solution [28, 33, 46]. The decoupling procedure enjoys the use of a local static condensation for the flow and mechanics. The MFE system resulting from each subsystem can be reduced to a symmetric and positive definite one; pressure is the sole unknown for the flow problem, and the displacement and rotation (may also be only the displacement depending on the used quadrature rule) for the mechanics (see [6, 7] for more details). Particularly, in the last years, a lot of research has been done on the fixed-stress method [5, 13, 18, 23]. Applied to problem (1.3), it can be rewritten, as, see [11]:

Definition 1.2** (The space-time fixed-stress algorithm).**

Chose an initial approximation ${\boldsymbol{\sigma}}^{0}\in H^{1}(0,T;\amsmathbb{W})$ of ${\boldsymbol{\sigma}}$ and a tolerance $\epsilon>0$. Set $k:=-1$. 2. 2.

Do**

1. (a)

Increase $k:=k+1$. 2. (b)

Compute $({\mathbf{w}}^{k},p^{k})\in L^{2}(0,T;\mathbf{W})\times H^{1}(0,T;Q)$ such that

[TABLE] 3. (c)

Compute $({\boldsymbol{\sigma}}^{k},{\mathbf{u}}^{k},{\boldsymbol{\zeta}}^{k})\in H^{1}(0,T;\amsmathbb{W})\times L^{2}(0,T;\mathbf{Q})\times L^{2}(0,T;\amsmathbb{Q}_{\textnormal{sk}})$ such that

[TABLE]

While* $\dfrac{\|({\boldsymbol{\sigma}}^{k},p^{k})-({\boldsymbol{\sigma}}^{k-1},p^{k-1})\|_{L^{2}(\Omega\times[0,T])}}{\|({\boldsymbol{\sigma}}^{k-1},p^{k-1})\|_{L^{2}(\Omega\times[0,T])}}\geq\epsilon$.*

The above method is the space-time fixed stress introduced first [13], in which we solve first the problem of flow over the whole space-time interval, then exchange the space-time information to solve the mechanics. This method is of interest for (i) the flexibility to use different time steps for flow and mechanics (ii) the advantage to derive error and a posteriori error analysis, permitting the use of adaptive asynchronous time-stepping (iii) the possibility to parallelize step 2.(c) of the algorithm. The classical fixed-stress algorithm in-space with four-field mixed formulation was analyzed in [46], where a priori convergence results are given. Other related works on this method can be found in [5, 14, 16, 23, 33] and the references therein.

In this paper, we are interested in designing adaptive versions of two common discretizations of the algorithm addressed in Definition 1.2, see Algorithm 3.1, 3.2 (standard), and Algorithm 4.2 (adaptive) below. To this aim, two iterative solution strategies for the Biot’s consolidation problem are presented; they are based on the above fixed stress iterative scheme, in which at each iteration, the space-time subsystems are solved sequentially using MFE methods in-space and with a backward Euler scheme in-time (cf. [46]). We constitute their adaptive counterpart upon the distinction of the different error components arising in the standard fixed stress algorithm, namely the spatial discretization, the temporal discretization, and the fixed stress iteration components. To arrive to this aim, we take ideas from [8, 22, 25, 29], for general a posteriori error techniques taking into account inexact iterative solvers, but most closely form [1], where a domain decomposition problem is solved via space-time iterative methods. Particularly, we will rely on [3, Theorem 6.2] where an energy-type-norm differences between the exact and the approximate pressure and displacement is shown to be bounded by the dual norm of the residuals. The developed adaptive fixed-stress algorithm is applicable on any locally conservative discretization for the two coupled subsystems, such as cell-centered finite volume scheme, multipoint mixed finite element, mimetic finite difference and hybrid high-order discontinuous Galerkin [12, 26, 37, 41]. It can also be extended to conforming methods using equilibrated flux and stress reconstructions ( cf. [40]).

In contrast to what is developed in [3], three additional features to be treated in this work; first is that the current setting targets inexact iterative coupling schemes for the Biot system and not monolithic solvers; second that the MFE methods here provides at each iteration of the coupling algorithms approximate flux and stress not balanced with the source terms; third is that the actual setting provides adaptive asynchronous time-stepping for the flow and mechanics problems. Here, we first show that the presented a posteriori error estimate delivers sharp bound (as reflected by moderate effectivity indices) for the actual energy-type error, and this at each iteration of the coupling algorithm. We also show how the overall error propagates between the flow and mechanics subproblems during the iterative process, and then to address the question of when to stop the iterations. This question was asked in [5, 11, 13], where the practitioners iterate between the two coupled subsystems until some fixed tolerance has been reached. The used stopping criterion is in fact mostly related to the algebraic error, i.e., the closeness of $({\boldsymbol{\sigma}}^{k},p^{k})$ to the convergent solution $({\boldsymbol{\sigma}}^{\infty},p^{\infty})$ is only taken into account without reference to the underlying continuous Biot’s problem (1.1). Here, by distinguishing the space, time and coupling errors, the adaptive stopping criterion for the iterative scheme that we propose instead is when the coupling error does not contribute significantly to the overall error. In grosso modo, the standard approach stops the iterations at some arbitrary tolerance, which hopefully is sufficiently accurate (but perhaps not!), while the approach based on error estimates stops the iterations at the correct time. Adaptive stopping criteria via a posteriori error estimates in the context of other model problems are treated in [1, 4, 19, 25], see also the references therein. Furthermore, the resulting algorithms involve tuning parameters that can be optimized (see [43]); the results show how a posteriori error estimates can help optimize these parameters. To the best of our knowledge, this combination of features in the adaptive fixed-stress algorithms is unique.

The paper is organized as follows. Section 2 fix the notation for temporal and spatial meshes and defines some relevant functional spaces. In Section 3, we present two common discretizations of Algorithm (1.4)-(1.5), by combining in-space two mixed finite elements for the flow and mechanical problems, and a backward Euler scheme in-time. As a posteriori error estimate has no meaning for piecewise constant functions, the MFE approximate pressure and displacement will be locally postprocessed in order to obtain improved approximations. In Section 4, we first introduce two major improvements to these two standards algorithms, by designing for each one, an adaptive stopping criterion, and a balancing criterion equilibrating the space and time error components using an adaptive asynchronous time-stepping. These enhancements are used to design adaptive versions of the fixed-stress schemes based on a posteriori error estimates. We then construct the needed ingredients for the a posteriori error estimates: Section 5, defines the $H^{1}(\Omega)$- and $\mathbf{H}({\textnormal{div}},\Omega)$ conforming reconstructions. In Section 6, these ingredients are used to bound an energy-type error in the pressure and displacement at each iteration of the coupling algorithm by a guaranteed and fully computable error estimate. This a posteriori estimate is then elaborated by distinguishing the fixed-stress iterative coupling error from the space and time error components. We also separate the pressure error components from those of displacement errors. We show numerical results in Section 7. Finally, a conclusion that highlights our developments is given in Section 8.

2 Notation

We introduce here the partition of $\Omega$, time discretization, notation, and function spaces; see [2] for a similar notation.

2.1 Partitions of the time interval $(0,T)$

The space-time iterative method we use supports asynchronous time grids for flow and mechanics. To this aim, the subscripts “f”, and “m” will be used throughout, to stand for flow and mechanics, respectively. For integer values $N_{\textnormal{f}}>0$, let $\left(\tau^{n}_{\textnormal{f}}\right)_{1\leq n\leq N_{\textnormal{f}}}$ denote a sequence of positive real numbers corresponding to the discrete flow time steps such that $T=\sum_{n=1}^{N_{\textnormal{f}}}\tau^{n}_{\textnormal{f}}$. Let $t^{0}_{\textnormal{f}}:=0$, and $t_{\textnormal{f}}^{n}:=\sum_{j=1}^{n}\tau^{j}_{\textnormal{f}},\ 1\leq n\leq N_{\textnormal{f}}$ be the discrete times for the flow problem. Let $I^{n}_{\textnormal{f}}:=(t^{n-1}_{\textnormal{f}},t^{n}_{\textnormal{f}}],\ 1\leq n\leq N_{\textnormal{f}}$. For the time stepping for the problem of mechanics, we will restrict ourselves to the case in which a fixed number of local flow time steps corresponds to one coarse mechanics time step. We suppose that $N_{\textnormal{f}}=\delta_{\textnormal{fm}}N_{\textnormal{m}}$, with $\delta_{\textnormal{fm}}$ and $N_{\textnormal{m}}$ are given positive integer values, where $\delta_{\textnormal{fm}}$ is the fixed number of local flow time steps within one coarse mechanics time step. We then let $\left(\tau^{\ell}_{\textnormal{m}}\right)_{1\leq\ell\leq N_{\textnormal{m}}}$ such that $T=\sum_{\ell=1}^{N_{\textnormal{m}}}\tau^{\ell}_{\textnormal{m}}$; for $1\leq\ell\leq N_{\textnormal{m}}$. We have then $\tau^{\ell}_{\textnormal{m}}=\sum_{n=(\ell-1)\delta_{\textnormal{fm}}+1}^{\ell\delta_{\textnormal{fm}}}\tau^{n}_{\textnormal{f}}$, and we let $t^{0}_{\textnormal{m}}:=0$, and $t_{\textnormal{m}}^{\ell}:=\sum_{j=1}^{\ell}\tau^{\ell}_{\textnormal{m}},\ 1\leq\ell\leq N_{\textnormal{m}}$ be the discrete times for the problem of mechanics; see [5] for a similar notation. We use $I^{\ell}_{\textnormal{m}}:=(t^{\ell-1}_{\textnormal{m}},t^{\ell}_{\textnormal{m}}],\ 1\leq\ell\leq N_{\textnormal{m}}$. For any sufficiently smooth function $v_{h\tau}$, we use the notation $v_{h}^{n}:=v_{h\tau}(\cdot,t^{n}_{\textnormal{f}})$, for all $0\leq n\leq N_{\textnormal{f}}$.

2.2 Partition of the domain $\Omega$

Let $\mathcal{T}_{h}$ be a simplicial mesh of $\Omega$, matching in the sense that for two distinct elements of $\mathcal{T}_{h}$ their intersection is either an empty set or their common vertex or edge. Let $h_{K}$ denote the diameter of $K\in\mathcal{T}_{h}$ and $h$ be the largest diameter of all triangle; $h:=\max_{K\in\mathcal{T}_{h}}h_{K}$. The set of vertices of the mesh $\mathcal{T}_{h}$ is denoted by $\mathcal{V}_{h}$, $\mathcal{V}_{h}^{\textnormal{int}}$ for the set of interior vertices, and $\mathcal{V}_{h}^{\textnormal{ext}}$ for the set of boundary vertices. For each ${\mathbf{a}}\in\mathcal{V}_{h}$, let $\mathcal{T}_{h}^{{\mathbf{a}}}$ denote the patch of the vertex ${\mathbf{a}}$, i.e., all the elements $K\in\mathcal{T}_{h}$ which share the vertex ${\mathbf{a}}$. We denote by $\omega_{{\mathbf{a}}}$ the corresponding open subset of $\Omega$.

2.3 Discrete function spaces

To approximate the flow subproblem (1.4), we let $Q_{h}\times\mathbf{W}_{h}\subset Q\times\mathbf{W}$ be the Raviart–Thomas–Nédélec mixed finite element spaces of order zero on the mesh $\mathcal{T}_{h}$ (cf. [9]):

[TABLE]

where ${\mathbf{RTN}}_{0}(K)$ denotes the lowest-order Raviart–Thomas–Nédélec finite-dimensional subspace associated with the element $K\in\mathcal{T}_{h}$.

To approximate the mechanics subproblem (1.5), we let $\mathbf{Q}_{h}\times\amsmathbb{W}_{h}\times\amsmathbb{Q}_{\textnormal{sk},h}\subset\mathbf{Q}\times\amsmathbb{W}\times\amsmathbb{Q}_{\textnormal{sk}}$ be the Arnold–Falk–Winther mixed finite elements with weakly symmetric stress for the lowest-order stresses on the mesh $\mathcal{T}_{h}$ (cf. [10]):

[TABLE]

where $[\amsmathbb{P}_{0}(K)]_{\textnormal{sk}}^{d\times d}$ denotes the subspace of $[\amsmathbb{P}_{0}(K)]^{d\times d}$ composed of skew symmetric–valued tensors.

Let $E$ be a space of functions defined on $\Omega$. We denote $P^{1}_{\tau}(E)$ the vector space of functions continuous in time and with values in $E$. We also denote by $P^{0}_{\tau}(E)$ the space of functions piecewise constant in time and with values in $E$. We have then if $v_{h\tau}\in P^{1}_{\tau}(E)$, then $\partial_{t}v_{h\tau}\in P^{0}_{\tau}(E)$ is such that for all $1\leq n\leq N_{\textnormal{f}}$,

[TABLE]

3 Fully discrete space-time fixed-stress schemes based on MFE in-space

and the backward Euler scheme in-time

In this section, we provide two discretization of Algorithm (1.4)-(1.5) using the backward Euler scheme in-time, and in-space, using two mixed finite elements methods for the linear elasticity and flow problems. A post-processing routine is then given to preview the numerical pressure and displacement solutions.

3.1 Two standard discrete fixed-stress schemes.

In the first algorithm, we consider the case of equal time grids for the flow and mechanics problems, i.e., $N_{\textnormal{f}}=N_{\textnormal{m}}$. The fully discrete form of Algorithm 1.2 reads then as follows:

Algorithm 3.1** (The global-in-time fixed-stress).**

Chose an initial approximation ${\boldsymbol{\sigma}}^{0}_{h\tau}\in P^{0}_{\tau}(\mathbf{W}_{h})$ of ${\boldsymbol{\sigma}}$, a real constant $\beta>0$, and a tolerance $\epsilon>0$. Set $k:=-1$. 2. 2.

Do**

1. (a)

Increase $k:=k+1$ and set $n:=0$. 2. (b)

Do**

1. i.

Increase $n:=n+1$. 2. ii.

Approximate $({\mathbf{w}}_{h}^{k,n},p^{k,n}_{h})\in\mathbf{W}_{h}\times Q_{h}$, the solution to

[TABLE]

While* $n\leq N_{\textnormal{f}}$.* 3. (c)

Reset $n:=0$. 4. (d)

Do**

1. i.

Increase $n:=n+1$. 2. ii.

Approximate $({\boldsymbol{\sigma}}_{h}^{k,n},{\mathbf{u}}_{h}^{k,n},{\boldsymbol{\zeta}}_{h}^{k,n})\in\amsmathbb{W}_{h}\times\mathbf{Q}_{h}\times\amsmathbb{Q}_{\textnormal{sk},h}$, solution to

[TABLE]

While* $n\leq N_{\textnormal{f}}$.*

While* $\left(\dfrac{\sum_{n=1}^{N_{\textnormal{f}}}\|({\boldsymbol{\sigma}}^{k,n}_{h},p^{k,n}_{h})-({\boldsymbol{\sigma}}^{k-1,n}_{h},p^{k-1,n}_{h})\|^{2}}{\sum_{n=1}^{N_{\textnormal{f}}}\|({\boldsymbol{\sigma}}^{k-1,n}_{h},p^{k-1,n}_{h})\|^{2}}\right)^{\frac{1}{2}}\geq\epsilon$. \tagform@3.3*

We present now the nonconforming-in-time counterpart of Algorithm 1.2 in the spirit of multi-rate fixed-stress scheme specified in [5]:

Algorithm 3.2** (The nonconforming-in-time (multi-rate) fixed-stress).**

Chose an initial approximation ${\boldsymbol{\sigma}}^{0}_{h\tau}\in P^{0}_{\tau}(\mathbf{W}_{h})$ of ${\boldsymbol{\sigma}}$, a real constant $\beta>0$, and a tolerance $\epsilon>0$. Set $\ell:=-\delta_{\textnormal{fm}}$. 2. 2.

Do**

1. (a)

Increase $\ell:=\ell+\delta_{\textnormal{fm}}$ and set $k:=-1$. 2. (b)

Do**

1. i.

Increase $k:=k+1$ and set $m:=0$. 2. ii.

Do**

1. A.

Increase $m:=m+1$. 2. B.

Approximate $({\mathbf{w}}_{h}^{k,\ell+m},p^{k,\ell+m}_{h})\in\mathbf{W}_{h}\times Q_{h}$, solution to

[TABLE]

While* $m\leq\delta_{\textnormal{fm}}$.* 3. iii.

Approximate $({\boldsymbol{\sigma}}_{h}^{k,\ell+\delta_{\textnormal{fm}}},{\mathbf{u}}_{h}^{k,\ell+\delta_{\textnormal{fm}}},{\boldsymbol{\zeta}}_{h}^{k,\ell+\delta_{\textnormal{fm}}})\in\amsmathbb{W}_{h}\times\mathbf{Q}_{h}\times\amsmathbb{Q}_{\textnormal{sk},h}$, solution to

[TABLE]

While* $\left(\dfrac{\|{\boldsymbol{\sigma}}^{k,\ell+\delta_{\textnormal{fm}}}_{h}-{\boldsymbol{\sigma}}^{k-1,\ell+\delta_{\textnormal{fm}}}_{h}\|^{2}+\sum_{m=1}^{\delta_{\textnormal{fm}}}\|p^{k,\ell+m}_{h}-p^{k-1,\ell+m}_{h}\|^{2}}{\|{\boldsymbol{\sigma}}^{k-1,\ell+\delta_{\textnormal{fm}}}_{h}\|^{2}+\sum_{m=1}^{\delta_{\textnormal{fm}}}\|p^{k-1,\ell+m}_{h}\|^{2}}\right)^{\frac{1}{2}}\geq\epsilon$. \tagform@3.6*

While* $\ell<\delta_{\textnormal{fm}}N_{\textnormal{m}}$.*

Remark 3.3** (The multi-rate FS).**

The convergence of the multi-rate fixed-stress was shown in [5] where mixed finite element method is used for the flow equations and where the mechanics is solved by conformal Galerkin method. Therein, the algorithm is also limited to one coarser time step for the mechanics and there is no study on the propagation of error due to temporal and spatial discretizations.

Remark 3.4** (Space-time vs multi-rate).**

We first notice that Algorithm 3.2 is practical to problems with a long time integration interval; In contrast to Algorithm 3.1, it requires reasonable computation ability and less storage resources to handle large-scale applications. Furthermore, Algorithm 3.2 is exploiting the different time scales for the flow and mechanics subsystems. We note also that the efficiency of the two algorithms can be improved when the free parameter $\beta$ is well-chosen (see [43]) and that step 2.(d) of Algorithm 3.1 in practice, is done in parallel as in [13].

3.2 Post-processing

We do here some improvements to the approximate solution $(p^{k}_{h\tau},{\mathbf{u}}^{k}_{h\tau})$ (cf. [21, 30]). This step is also mandatory to design from Algorithm 3.1 and 3.2, their adaptive versions based on energy-norm-type a posteriori error estimate. This is customary in mixed finite elements schemes, as an energy-norm-type a posteriori error estimate has no meaning for the piecewise constant, i.e., $\nabla p^{k,n}_{h}=\nabla{\mathbf{u}}^{k,n}_{h}=0$.

Let us notice first that in Algorithm 3.2, the approximate solution $({\mathbf{u}}^{k}_{h\tau},\boldsymbol{\sigma}^{k}_{h\tau})$ of the mechanics problem is defined in different time grids from the approximate flow solution $(p^{k}_{h\tau},{\mathbf{w}}^{k}_{h\tau})$, so we cannot proceed to the post-processing of the displacement and the reconstruction of the stress tensor unless we build the couple $({\mathbf{u}}^{k}_{h\tau},\boldsymbol{\sigma}^{k}_{h\tau})$ at the finer time steps $t^{n}_{\textnormal{f}}=t^{\ell+m}_{\textnormal{f}}$, for all $\ell=0,\delta_{\textnormal{fm}},2\delta_{\textnormal{fm}},3\delta_{\textnormal{fm}},\cdots,(N_{\textnormal{m}}-1)\delta_{\textnormal{fm}}$, and for all $1\leq m\leq\delta_{\textnormal{fm}}-1$. To this aim, we construct the displacement and the stress tensor as follows: for $\ell=0,\delta_{\textnormal{fm}},2\delta_{\textnormal{fm}},3\delta_{\textnormal{fm}},\cdots,(N_{\textnormal{m}}-1)\delta_{\textnormal{fm}}$, we set

[TABLE]

Note that this post-processing is explicit and its cost is negligible. The post-processing of the pressure $p^{k,n}_{h}$ is as follows [21, 30, 24]: at each iteration $k\geq 1$, we calculate the improved solution $\widetilde{p}^{k,n}_{h}\in\amsmathbb{P}_{2}(\mathcal{T}_{h})$ in each element $K\in\mathcal{T}_{h}$ such that

[TABLE]

This post-processing is computationally cheap and easy to be implemented. We then extend this post-processing (cf. [30]) to the vector case, leading to define a function $\widetilde{{\mathbf{u}}}^{k,n}_{h}\in\left[\amsmathbb{P}_{2}(\mathcal{T}_{h})\right]^{d}$, such that

[TABLE]

where $\mathbf{e}_{i}\in{\amsmathbb R}^{d}$ denotes the $i$-th Euclidean unit vector. We finally assume for simplicity that the initial conditions are satisfied exactly, i.e.,

[TABLE]

We define the continuous and piecewise affine in–time functions $\widetilde{p}^{k}_{h\tau}$ and $\widetilde{{\mathbf{u}}}^{k}_{h\tau}$ by

[TABLE]

The key observation in the above post-processing is that they use only local operations, which are independent from each other and hence parallelizable.

4 The adaptive fixed-stress algorithms

The purpose of this section is to reduce as much as possible the computational effort of Algorithm 3.1 and 3.2 as in [1, 4, 20, 22]. The improvements of these two standards algorithms stems from (i) important savings in terms of the number of coupling iterations can be achieved using adaptive stopping criterion (ii) a significant gain in the computational resources is obtained by balancing the error components via an asynchronous adaptivity of the temporal meshes (iii) optimizing the tuning parameter $\beta$.

4.1 Methodology for adaptive asynchronous time-stepping and adaptive stopping criteria

Let $\eta^{k,n}_{\textnormal{sp,J}}$, $\eta^{k,n}_{\textnormal{tm,J}}$ and $\eta^{k,n}_{\textnormal{it,J}}$, for $\textnormal{J}=\textnormal{P,\,U}$, be respectively the estimators of the spatial discretization error, the temporal discretization error and the fixed-stress coupling error at the $n\textnormal{-th}$ time step and on the $k\textnormal{-th}$ iteration, where the index J=P is for the pressure error components, and that of J=U is for the displacement error components. We let $\eta^{k,n}_{\bullet}:=\eta^{k,n}_{\bullet,\textnormal{P}}+\eta^{k,n}_{\bullet,\textnormal{U}}$.

The first step of our developments is to equip Algorithm 3.1 and 3.2 with adaptive asynchronous time-stepping. To this aim, we propose to equilibrate the time errors with the spatial errors as follows; we adjust the time steps $\tau^{n}_{\textnormal{f}}$ and $\tau^{n}_{\textnormal{m}}$ so that

[TABLE]

where $\gamma_{\textnormal{tm,J}}$ and $\Gamma_{\textnormal{tm,J}}$, $\textnormal{J=P,\,U},$ are user-given weights, typically close to $1$.

An alternative to (4.1) being to balance the time errors from the flow and mechanics discretization with the global error by selecting the time steps $\tau^{n}_{\textnormal{m}}$ and $\tau^{n}_{\textnormal{f}}$ in such a way that

[TABLE]

The balancing criterion (4.1) controls the contributions of $\eta^{k,n}_{\textnormal{tm,P}}$ and $\eta^{k,n}_{\textnormal{tm,U}}$ in the overall error and leads to $\eta^{k,n}_{\textnormal{sp}}\approx\eta^{k,n}_{\textnormal{tm}}$. That of the second criterion (4.2) leads to equilibrate the time errors from the flow and mechanics, i.e., $\eta^{k,n}_{\textnormal{tm,P}}\approx\eta^{k,n}_{\textnormal{tm,U}}$.

The second step of our developments is to equip Algorithm 3.1 and 3.2 with adaptive stopping criteria. We then introduce a real parameter $\gamma_{\textnormal{it}}$ to be given in $(0,1)$. The stopping criteria for Algorithm 3.1 is chosen at each iteration $k$ as

[TABLE]

which implies that we stop the iterations if the coupling error is sufficiently lower than one of the other components. The stopping criteria for Algorithm 3.2 is set similarly: at each coarse mechanics time step $\ell=0,\delta_{\textnormal{fm}},2\delta_{\textnormal{fm}},\cdots,(N_{\textnormal{m}}-1)\delta_{\textnormal{fm}}$, at each iteration $k\geq 1$,

[TABLE]

Remark 4.1** (Algebraic errors).**

The systems within the flow and mechanics subsystems are solved with direct solvers. The present adaptive approach can also be combined with an iterative solver for each subproblem, and to further save computational effort, these latter can be stopped whenever the algebraic errors does not contribute significantly to the overall error, following [25].

4.2 The adaptive algorithms

We are now ready to present the adaptive counterparts of Algorithm 3.1 and 3.2, i.e., they are now equipped with adaptive asynchronous time-stepping and a posteriori stopping criterion. The adaptive version of Algorithm 3.1 is as follows:

Algorithm 4.2** (Adaptive Fixed-Stress with Asynchronous Time Mesh Refinement and Adaptive Stopping Criteria ).**

In step 1, chose also a real parameter $\gamma_{\textnormal{it}}$ and the real weights $\gamma_{\textnormal{tm}}$ and $\Gamma_{\textnormal{tm}}$, set $t_{\textnormal{f}}^{0}=0$, and give an initial ratio $\delta_{\textnormal{fm}}^{0}$ and time step for the flow $\tau^{0}_{\textnormal{f}}$ and the temporal refinement threshold $\tau_{\textnormal{min}}$. Set $k:=-1$. 2. 2.

Do**

1. (a)

Increase $k:=k+1$ and set $n:=0$. 2. (b)

Do**

1. i.

Increase $n:=n+1$. 2. ii.

Set $\tau^{n}_{\textnormal{f}}:=\tau^{n-1}_{\textnormal{f}}$. 3. iii.

Approximate $({\mathbf{w}}_{h}^{k,n},p^{k,n}_{h})$ by (3.1). 4. iv.

Calculate the estimators $\eta^{k,n}_{\textnormal{sp,P}}$, $\eta^{k,n}_{\textnormal{tm,P}}$ and $\eta^{k,n}_{\textnormal{it,P}}$. Check the balancing criterion (4.1) (or (4.2)). If not satisfied, refine or redefine the flow time step $\tau^{n}_{\textnormal{f}}$ in such a way that condition (4.1) (or (4.2)) holds or $\tau^{n}_{\textnormal{f}}\leq\tau_{\textnormal{min}}$, and return to step 2.(b)iii. 5. v.

Set $t^{n}_{\textnormal{f}}:=t^{n-1}_{\textnormal{f}}+\tau^{n}_{\textnormal{f}}$.

While* $t_{\textnormal{f}}^{n}\leq T$.* 3. (c)

Reset $n:=0$ and $t_{\textnormal{m}}^{0}=0$, and let $\tau^{0}_{\textnormal{m}}=\delta_{\textnormal{fm}}^{0}\tau^{0}_{\textnormal{f}}$. 4. (d)

Do**

1. i.

Increase $n:=n+1$. 2. ii.

Set $\tau^{n}_{\textnormal{m}}:=\tau^{n-1}_{\textnormal{m}}$. 3. iii.

Approximate $({\boldsymbol{\sigma}}_{h}^{k,n},{\mathbf{u}}_{h}^{k,n},{\boldsymbol{\zeta}}_{h}^{k,n})$ by (3.2). 4. iv.

Calculate the estimators $\eta^{k,n}_{\textnormal{sp,U}}$, $\eta^{k,n}_{\textnormal{tm,U}}$ and $\eta^{k,n}_{\textnormal{it,U}}$. Check the balancing criterion (4.1) (or (4.2)). If not satisfied, redefine the time step for the mechanics $\tau^{n}_{\textnormal{m}}$ using using $\left(\tau^{n}_{\textnormal{f}}\right)_{n\geq 0}$ in such a way that condition (4.1) (or (4.2)) holds, and return to step 2.(d)iii. 5. v.

Set $t^{n}_{\textnormal{m}}:=t^{n-1}_{\textnormal{m}}+\tau^{n}_{\textnormal{m}}$.

While* $t_{\textnormal{m}}^{n}\leq T$.*

Until* the criteria (4.3) is satisfied.*

Similarly, we propose to modify Algorithm 3.2. This yields to an adaptive fixed-stress scheme equivalent to Algorithm 4.2 but applied through temporal windowing technique. Precisely, the whole time interval $[0,T]$ is now split into $N$ time-windows $[0,T_{1}]$, $[T_{1},T_{2}],\cdots,[T_{N-1},T]$. Algorithm 4.2 is first applied on the first time window $[0,T_{1}]$. Afterwards, one applies the algorithm on the next time window $[T_{i-1},T_{i}]$ imposing as initial condition for $t^{0}_{\textnormal{f}}=t^{0}_{\textnormal{m}}=T_{i-1}$ the solution of the converged iterate of the end of the previous time window and proceeds in such a way until all time windows have been treated.

Remark 4.3** (Space adaptivity).**

The estimators are calculated on each element of the mesh and on each time step, and could also be used as indicators in order to refine adaptively the spatial mesh $\mathcal{T}_{h}$, so that the local spatial error estimators are distributed equally; see [3, 20, 44] and the references therein. Furthemore, the efficiency of the adaptive algorithms can be enhanced by using adaptive multiscale meshes, where two spatial meshes for the flow and mechanics subsystems are considered and where they are refined/coarsened adaptively in order to equilibrate the space errors for the two subsystems; see the multiscale discretizations techniques in [18, 38].

Remark 4.4** (Static condensation).**

The decoupling procedure permits the use of a local static condensation for the flow and mechanics and then to reduce the MFE system resulting from each subproblem to a symmetric and positive definite one; for the pressure for the flow problem, and for the displacement and rotation for the mechanics problem; with the same way as in [6, 7]. These systems are smaller and easier to solve than the original saddle point problems, but no further reduction is possible.

Remark 4.5** (Computational cost).**

In practice, the adaptive time-stepping strategy is only done in the first or second iterations of the Algorithm. Further, when the step size is modified once at current time step, the updated step size can be used for some subsequent time steps, say for example $5$ time steps. Then on the sixth time step, step 2.(b)iv is checked again if the step size needs to be modified for the next $5$ time steps. Furthermore, steps 2.(b)iv and 2.(d)iv can be done only every few iterations of the Algorithm.

5 Concept of $H^{1}$-, and $\mathbf{H}(\textnormal{div})$-conforming reconstructions in Biot’s poro-elasticity system

In this section, we develop basic tools that will allow us to build the estimators involved in the adaptive fixed-stress algorithms.

5.1 Pressure and displacement reconstructions

We first construct from $\widetilde{p}^{k,n}_{h}$, a $H^{1}_{0}$-conforming function $\hat{p}^{k,n}_{h}$ satisfying the mean value constraint

[TABLE]

To this aim, we proceed as in [21]; from the available post-processed pressure $\widetilde{p}^{k,n}_{h}$ at each iteration $k\geq 1$, we set

[TABLE]

Here ${\mathbf{a}}$ are the Lagrangian nodes situated in the interior of $\Omega$, $b_{K}$ denotes the standard (time-independent) bubble function supported on $K$, defined as the product of the barycentric coordinates of $K$, for all $K\in\mathcal{T}_{h}$, and scaled so that its maximal value is 1, and $\mathcal{I}_{\textnormal{av}}:\amsmathbb{P}_{2}(\mathcal{T}_{h})\rightarrow\amsmathbb{P}_{2}(\mathcal{T}_{h})\cap H^{1}(\Omega)$ is the interpolation operator given by

[TABLE]

At the Lagrange nodes ${\mathbf{a}}$ situated on the boundary $\partial\Omega$, we set ${\hat{p}}^{k,n}_{h}({\mathbf{a}}):=0$. In order to guarantee that the mean value constraint (5.1) holds true, we choose

[TABLE]

The same procedure can be applied to the post-processed displacement $\widetilde{{\mathbf{u}}}^{k,n}_{h}\in\left[\amsmathbb{P}_{2}(\mathcal{T}_{h})\right]^{d}$. This leads to a $\mathbf{H}^{1}_{0}(\Omega)$-conforming vector function $\hat{{\mathbf{u}}}^{k,n}_{h}$, satisfying the following mean value constraint:

[TABLE]

We end up with continuous and piecewise affine in–time functions $\widetilde{p}^{k}_{h\tau}$ and $\widetilde{{\mathbf{u}}}^{k}_{h\tau}$ by setting

[TABLE]

An interesting result of the above reconstructions is given in the following Lemma (cf. [3]):

Lemma 5.1** (Properties of $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})$).**

At each iteration $k\geq 1$ of Algorithm 4.2, let $(\widetilde{p}^{k}_{h\tau},\widetilde{{\mathbf{u}}}^{k}_{h\tau})$ be the post-processed pressure and displacement, and $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})$ be the reconstructed pressure and displacement. Then, for all $1\leq n\leq N_{\textnormal{f}}$, there holds

[TABLE]

5.2 Equilibrated flux $\hat{{\mathbf{w}}}^{k}_{h\tau}$ and

stress $\hat{\boldsymbol{\sigma}}^{k}_{h\tau}$ reconstructions

The second ingredient for the derivation of our a posteriori error estimates is to reconstruct an equilibrated flux $\hat{{\mathbf{w}}}^{k}_{h\tau}\in P^{0}_{\tau}(\mathbf{H}({\textnormal{div}},\Omega))$, locally conservative on the mesh $\mathcal{T}_{h}$:

[TABLE]

and to reconstruct an equilibrated stress $\hat{\boldsymbol{\sigma}}^{k}_{h\tau}\in P^{0}_{\tau}(\amsmathbb{H}({\textnormal{div}},\Omega))$, locally conservative on the mesh $\mathcal{T}_{h}$:

[TABLE]

These reconstructions are based on solving local Neumann problems by mixed finite elements posed over patches of elements around mesh vertices (cf. [21, 22, 40]). For each vertex ${\mathbf{a}}\in\mathcal{V}_{h}$, we introduce the mixed Raviart–Thomas finite element spaces posed on the patch domain $\omega_{{\mathbf{a}}}$:

[TABLE]

We then introduce the Arnold–Falk–Winther mixed finite elements spaces, posed on the patch domain $\omega_{{\mathbf{a}}}$, for all ${\mathbf{a}}\in\mathcal{V}_{h}$:

[TABLE]

We obtain the equilibrated velocity field $\hat{{\mathbf{w}}}^{k}_{h\tau}$, by solving first for $(\hat{{\mathbf{w}}}_{{\mathbf{a}}}^{k,n},q_{{\mathbf{a}}}^{k,n})\in\mathbf{W}^{{\mathbf{a}}}_{h}\times Q^{{\mathbf{a}}}_{h}$, for all $1\leq n\leq N_{\textnormal{f}}$, such that

[TABLE]

For the equilibrated stress $\hat{\boldsymbol{\sigma}}_{h\tau}^{k}$, we solve local Neumann mechanics problems by mixed finite elements, with weakly symmetric stress: find $(\hat{\boldsymbol{\sigma}}_{{\mathbf{a}}}^{k,n},{\mathbf{z}}_{{\mathbf{a}}}^{k,n},\boldsymbol{\beta}_{{\mathbf{a}}}^{k,n})\in\amsmathbb{W}^{{\mathbf{a}}}_{h}\times\mathbf{Q}^{{\mathbf{a}}}_{h}\times\amsmathbb{Q}_{\textnormal{sk},h}^{{\mathbf{a}}}$, for all $1\leq n\leq N_{\textnormal{f}}$, such that

[TABLE]

The above local problems are well-posed owing to the properties of mixed finite elements (cf. [9, 10]). We can easily observe that the computational cost of the flux and stress reconstructions can be substantially reduced by pre-processing, a step that is fully parallelizable.

6 The a posteriori error estimates

We derive in this section, based on the previous constructions, a posteriori error estimates for the solution of Algorithm 3.1 or 3.2. This is done by bounding an energy error between the exact weak solution $(p,{\mathbf{u}})$ of problem (1.1) and the approximate solution $(\widetilde{p}^{k}_{h\tau},\widetilde{{\mathbf{u}}}^{k}_{h\tau})$, by a guaranteed and fully computable upper bound, and this at each iteration $k\geq 1$ of Algorithm 3.1 or 3.2.

6.1 The error measure

The first question in a posteriori error estimates is that of the error measure; here we will in particular rely on  [3, Theorem 6.2], where an energy-type error in the pressure and displacement is shown to be bounded by the dual norm of the residuals, and where the Biot’s consolidation equations (1.3), (1.1c) are discretized using MFE method in-space and with a backward Euler scheme in-time and solved monolithically. To this aim, and for all times $t\in(0,T]$, we let

[TABLE]

and introduce the energy space

[TABLE]

Then, we introduce a weak formulation of (1.1): find $(p,{\mathbf{u}})\in\mathcal{E}_{T}$ such that $p(\cdot,0)=p_{0}$ and ${\mathbf{u}}(\cdot,0)={\mathbf{u}}_{0}$ and such that

[TABLE]

where $\langle\cdot,\cdot\rangle$ denotes the duality pairing between $H^{-1}(\Omega)$ and $H^{1}_{0}(\Omega)$. The existence and uniqueness of the solution to this problem was addressed in [3]. Still following [3], we introduce the following energy-type error measure

[TABLE]

The above norms are well-defined owing to the properties of the weak solution $(p,{\mathbf{u}})$ and the reconstructed functions $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})$, i.e., we have both $(p,{\mathbf{u}})$ and $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})$ in $\mathcal{E}_{T}$.

6.2 The error estimators

Before formulating our estimators, we define the broken Sobolev space $H^{1}(\mathcal{T}_{h})$ as the space of all functions $v\in L^{2}(\Omega)$ such that $v|_{K}\in H^{1}(K)$, for all $K\in\mathcal{T}_{h}$. The energy semi-norm on $H^{1}(\mathcal{T}_{h})$ is given by

[TABLE]

where the sign $\nabla$ denote the element-wise gradient, i.e., the gradient of a function restricted to a mesh element $K\in\mathcal{T}_{h}$. The energy norm in $\mathbf{L}^{2}(\Omega)$ is given by

[TABLE]

We also recall the Poincaré inequality:

[TABLE]

where $q_{K}$ is the mean value of the function $q$ on the element $K$ given by $q_{K}=\int_{K}q\,\textnormal{d}{\mathbf{x}}/|K|$ and $C_{P,K}=1/\pi$ whenever the element $K$ is convex. In what follows, we denote respectively by $c_{\mathbf{K},K}$ and $C_{\mathbf{K},K}$ the smallest and the largest eigenvalue of the tensor $\mathbf{K}$ in $K\in\mathcal{T}_{h}$. We introduce the local residual estimators

[TABLE]

the flux estimators

[TABLE]

the pressure nonconformity estimators

[TABLE]

Therefrom, we introduce the global versions by

[TABLE]

6.3 Guaranteed reliability

We now provide a guaranteed estimate on the total error in particular valid on each iteration of Algorithm 3.1 or 3.2. This result extend the results of our previous work [3, Section 5], where a computable guaranteed bound on the energy-type error between the exact solution and its approximation with an exact solver has been derived.

Theorem 6.1** (Global-in-time a posteriori error estimate).**

Let $(p,{\mathbf{u}})$ be the weak solution of (6.1). At an arbitrary iteration $k\geq 1$ of Algorithm 3.1 or 3.2, let $(\widetilde{p}^{k}_{h\tau},\widetilde{{\mathbf{u}}}^{k}_{h\tau})$ be the post-processed pressure and displacement of subsection 3.2, $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})$ be the reconstructed pressure and displacement and $(\hat{{\mathbf{w}}}^{k}_{h\tau},\hat{\boldsymbol{\sigma}}^{k}_{h\tau})$ be the reconstructed flux and stress of section 5. Then, there holds

[TABLE]

where

[TABLE]

Notice that we have set $L_{\textnormal{P}}=1$ and $L_{\textnormal{U}}=\frac{1}{\mu}$, and for $1\leq n,l\leq N_{\textnormal{f}}$,

[TABLE]

Proof.

Recalling (6.2a), and (3.10), we have from [3, Theorem 6.2], for any given couple $(\hat{p}_{h\tau},\hat{{\mathbf{w}}}_{h\tau})\in\mathcal{E}_{T}$,

[TABLE]

featuring the residuals $\mathcal{R}_{\textnormal{P}}(\hat{p}_{h\tau},\hat{{\mathbf{u}}}_{h\tau})\in X^{\prime}_{T}$ and $\mathcal{R}_{\textnormal{U}}(\hat{p}_{h\tau},\hat{{\mathbf{u}}}_{h\tau})\in\mathbf{X}^{\prime}_{T}$ of the weak formulation (6.1): for all $q\in X_{T}$ and ${\mathbf{v}}\in\mathbf{X}_{T}$,

[TABLE]

The dual norms of the residuals are given by

[TABLE]

At each iteration $k\geq 1$, the approximate solution $(\widetilde{p}^{k}_{h\tau},\widetilde{{\mathbf{u}}}^{k}_{h\tau})$ is not an element of $\mathcal{E}_{T}$, contrarily to the reconstructed solution, i.e., $(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})\in\mathcal{E}_{T}$. Thus, to use (6.12), we apply the triangle inequality to get

[TABLE]

where we can bound the first term of the right-hand side using (6.12). What remains is to give a computable upper bound for the residuals $||\mathcal{R}_{\textnormal{P}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||_{X^{\prime}_{T}}$ and $||\mathcal{R}_{\textnormal{U}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||_{\mathbf{X}^{\prime}_{T}}$ together with $||(\hat{p}^{k}_{h\tau}-\widetilde{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau}-\widetilde{{\mathbf{u}}}^{k}_{h\tau})||_{\textnormal{en}}$ and then combine these results.

1) A computable upper bound for $||\mathcal{R}_{\textnormal{P}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||_{X^{\prime}_{T}}$ and $||\mathcal{R}_{\textnormal{U}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||_{\mathbf{X}^{\prime}_{T}}$. Proceeding as in [2, 40], adding $(\hat{{\mathbf{w}}}^{k,n}_{h},\nabla q)$ to (6.13), choosing $q\in X_{T}$ with $||q||_{X_{T}}=1$ and applying the Green theorem, and using (5.6), we obtain

[TABLE]

Then, it can be inferred using (6.15) and the Poincaré inequality (6.5) followed by Cauchy-Schwarz inequality that

[TABLE]

We repeat the same steps for (6.14), by adding and subtracting $(\hat{{\boldsymbol{\sigma}}}^{k,n}_{h\tau},\nabla{\mathbf{v}})$ (we replace $(\hat{{\boldsymbol{\sigma}}}^{k,n}_{h\tau},\boldsymbol{\epsilon}({\mathbf{v}}))$ by $(\hat{{\boldsymbol{\sigma}}}^{k,n}_{h\tau},\nabla{\mathbf{v}})$ due to symmetry), using (5.7) and applying the Poincaré inequality (6.5) together with the Cauchy-Schwarz inequality and definition (6.16), we obtain

[TABLE]

Replacing (6.19) and (6.20) in (6.12), so we are left to bound the third and fourth terms of the right-hand side of 6.12. Using the fact that $||\mathcal{R}_{\textnormal{U}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||_{\mathbf{X}_{t}^{\prime}}^{2}$ (also $||\mathcal{R}_{\textnormal{P}}(\hat{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau})||^{2}_{X^{\prime}_{t}}$) is a nondecreasing function of the time $t$, we easily obtain from (6.19)-(6.20),

[TABLE]

In a similar way, we infer

[TABLE]

and similarly

[TABLE]

We use (6.19)–(6.23) in (6.12), thus we bound the first term of (6.17):

[TABLE]

2) A computable upper bound to $||(\hat{p}^{k}_{h\tau}-\widetilde{p}^{k}_{h\tau},\hat{{\mathbf{u}}}^{k}_{h\tau}-\widetilde{{\mathbf{u}}}^{k}_{h\tau})||_{\textnormal{en}}$. To bound this term presenting the nonconformity estimator, we proceed as in  [3, Theorem 5.3 $\&$ Lemma 5.7], we promptly arrive to

[TABLE]

The estimate (6.10) is obtained by replacing (6.24) and (6.25) in (6.17). ∎

6.4 An a posteriori error estimate distinguishing the space, time and fixed-stress coupling

errors

Our goal in this section is to distinguish the different error components. Particularly, we separate the iterative coupling error from the estimated space and time errors, which are predefined and used efficiently in Algorithm 4.2. To this purpose, at each iteration $k\geq 1$, we define for all $1\leq n\leq N_{\textnormal{f}}$, the local spatial, temporal and iterative coupling estimators

[TABLE]

Therefrom, we introduce like in (6.9), $\textnormal{a}=\textnormal{sp,\,tm,\,it}$,

[TABLE]

and then introduce their global versions like in (6.11) by

[TABLE]

where $\delta_{\textnormal{a}}=0$ for $\textnormal{a}=\textnormal{tm,\,it}$ and $\delta_{\textnormal{sp}}=1$.

Remark 6.2** (Fixed-stress estimator).**

In the above estimators, we have to mention that, for conforming time discretization, the iterative coupling estimators $\eta^{k}_{\textnormal{it,J}}$, J=U, P, tends to zero when the fixed-stress algorithm converges. However, this is not true for non-conforming time discretization as in the multi-rate algorithm or the adaptive one. Precisely, at each iteration $k\geq 1$, the reconstructed flux and stress are satisfying respectively the local mass conservation (5.10) and (5.11) which is not the case for the approximate flux from (3.4) and the approximate stress from (3.5) (step 2.(b)ii in Algorithm 3.2 and step 2.(b)iii in Algorithm 4.2) and this is even when the fixed-stress converges. In other words, the fixed-stress estimator $\eta^{k}_{\textnormal{it}}=\eta^{k}_{\textnormal{it,P}}+\eta^{k}_{\textnormal{it,U}}$ becomes a non-conformity in-time estimator when the fixed stress algorithm converges.

Lemma 6.3** (A posteriori error estimate distinguishing error components).**

Let the assumptions of Theorem 6.1 be satisfied. Then there holds

[TABLE]

Proof.

We add and subtract the discrete fluxes in the flux estimator (6.7a) then apply the triangle inequality yield, for all $K\in\mathcal{T}_{h},\,\,t\in I_{\textnormal{f}}^{n}$,

[TABLE]

Another triangle inequality leads to

[TABLE]

In the same way, we obtain for the stress estimator (6.7b), for all $K\in\mathcal{T}_{h},\,\,t\in I_{\textnormal{f}}^{n}$,

[TABLE]

In these two inequalities, the first terms in the right-hand side form the contributions of the pressure and displacement in the fixed-stress error, whereas the second ones contribute to the two components of the space discretization error. The last terms can be integrated in time, yielding to the two components of the time error. What remains is to replace (6.31) and (6.32) in (6.10) with the use of (6.26)–(6.28), where we used the equality of norms in (6.31), i.e., $||\mathbf{K}\nabla v||_{\star,K}=||\mathbf{K}^{-\frac{1}{2}}(\mathbf{K}\nabla v)||_{K}=|||v|||_{K}$, leading to the estimate (6.29). ∎

7 Numerical results

In this section we illustrate the efficiency of our theoretical results on numerical experiments. We have chosen two examples designed to show how the adaptive fixed-stress scheme behaves vs the standard ones and this is done on different physical and geometrical situations.

7.1 Test problem 1: an academic example

with a manufactured solution

We consider in the computational domain $\Omega=[0,1]^{2}$ and the final time $T=1$. The analytical solution of Biot’s consolidation problem is supposed to be:

[TABLE]

The effective parameters are $\mathbf{K}={\mathbf{I}}$, and $\alpha=\mu=\lambda=1$. This analytical solution which has homogeneous initial and Dirichlet boundary values for $p$ and ${\mathbf{u}}$, generates from (1.1a)- (1.1b) non-zeros source terms $\mathbf{f}({\mathbf{x}},t)$ and $g({\mathbf{x}},t)$.

7.1.1 Stopping criteria balancing the error components

The aim here is to illustrate the performance of the adaptive stopping criteria introduced in Section 4.1. To this purpose, we consider a uniform space-time mesh with $h=1/16$, and $\tau_{\textnormal{f}}^{n}=\tau_{\textnormal{m}}^{n}=(2h)^{2}$. The tuning parameter is chosen $\beta=\dfrac{\alpha^{2}}{\delta(\frac{2\mu}{d}+\lambda)}$, with $\delta=2$. The choice of the parameter $\delta$ is theoretically calculated in [11, 43, 32] and possibly should lead to the best performance of the fixed-stress method in terms of number of iterations. We first test the performance of the space–time fixed stress algorithm (Algo. 3.1) equipped with an adaptive stopping criteria. Therein, we set $\gamma_{\textnormal{it}}=0.2$ and compare the results with the standard approach in which the fixed-stress algorithm is continued until the algebraic residual-based criteria (2) is satisfied for an (arbitrary) threshold $\epsilon$.

Figure 1(a) displays the dependence of the total error and of the various estimators on the fixed-stress iterations, where various stopping criteria are used. We can see that the space (blue) and time (green) estimators remain constant during the computation in contrast to the fixed-stress estimator, which gives a numerical indication that the error components are distinguished. For the fixed-stress estimator (red), we can see that the adaptive stopping criteria (4.3) is satisfied after 16 iterations only, while the classical one (2) needs 34 iterations with $\epsilon=10^{-6}$, 18 iterations with $\epsilon=10^{-3}$, and only 10 iterations with $\epsilon=10^{-2}$. We can remark that the total error (magenta) and the total estimator (black) decrease rapidly for the first 12 fixed-stress iterations, after which they decrease very slowly, as the influence of the fixed-stress iteration error becomes negligible. This is exactly the point where our adaptive fixed-stress method stops. This results in a significant saving of fixed-stress iterations as well as excludes possible inaccurate solutions from the algorithm (like with $\epsilon=10^{-2})$.) As an example, we make a gain of $53\%$ of total fixed-stress iterations compared to the classical algorithm with $\epsilon=10^{-6}$.

As known, the number of iterations to achieve convergence for the fixed-stress can differ considerably depending on the choice of the tuning parameter $\delta$. Thus, in Figure 1(b), we plot the number of iterations required by the fixed-stress algorithm as function of the parameter $\delta$. Therein, we stop the algorithm when the adaptive stopping criteria is satisfied. Clearly, the estimator behaves very similarly to what is usually observed for the fixed-stress error (see, e.g., [43]). Moreover, the theoretical parameter (marked by a star) coincides with the numerically optimal value. In Figure 2, we display the spatial distribution of the different error components (left) and of the corresponding estimators (right) at the final time $t=T$. Clearly, the distribution of the estimated errors reflects the exact ones. Also as expected, we observe in Figure 2(d) that the estimated fixed-stress error is sufficiently small to not contribute significantly in the overall error.

In Figure 3, the effectivity index for the space–time fixed stress approach is presented. It is calculated by the ratio between the total estimator and the exact total error at the iteration $k$ of the fixed-stress algorithm. The effectivity index oscillates during the first 10 iterations, then decreases to approximately 4.27 and then remains constant until the end of the computation. One of the reasons of why this factor is far from 1, may be that the negative norms involved in the exact error $\|(p-\tilde{p}^{k}_{h\tau},{\mathbf{u}}-\tilde{{\mathbf{u}}}^{k}_{h\tau})\|_{\textnormal{en}}$ are not calculated. Another explanation is that in practice, we use estimate (6.29) instead of (6.10) where the different error components are not yet separated.

In the second set of experiments, we study the performance of the adaptive stopping criteria on the multi-rate fixed-stress algorithm (MFS). Here, we compare the results with the classical multi-rate algorithm in which the algebraic residual-based criteria (2b) is used with various threshold $\epsilon$.

In Figures 4(a)–4(c), we plot the evolution of the total error and the various estimators on the fixed-stress iterations for the final coarse mechanics step. There, we compare the adaptive to the standard multi-rate algorithm for various ratios of discretization in time, $\delta_{\textnormal{fm}}=8,\,4,\,2$. We remark that (i) the discretization in space estimator (blue) as well as the discretization in time estimator (green) are approximately constant in each case (ii) the discretization in time estimator goes up and approaching from the discretization in space estimator when we increase the ratio $\delta_{\textnormal{fm}}$. Note that the discretization in space estimator is the same for the different ratios if it is scaled with the time steps. These findings confirm numerically that we have practically distinguished the time discretization error from the spatial discretization error. Concerning the fixed-stress estimator, we recall that in that case, $\eta^{k,\ell}_{\textnormal{it}}$ mixes fixed-stress and nonconformity-in-time errors (see Remark 6.2). Thus, we observe for the case $\delta_{\textnormal{fm}}=4$, and $\delta_{\textnormal{fm}}=8$ that, $\eta^{k,\ell}_{\textnormal{it}}$ dominates the total error until iteration 2 or 3, then becomes smaller than $\eta^{k,\ell}_{\textnormal{sp}}$ and $\eta^{k,\ell}_{\textnormal{tm}}$ until iteration 4 or 5, and therefrom remains constant as the influence of the fixed-stress iteration error becomes negligible compared to the nonconformity-in-time error. For $\delta_{\textnormal{fm}}=2$, the nonconformity error is small enough so as not to contribute in $\eta^{k,\ell}_{\textnormal{it}}$ until convergence. For the cases, $\delta_{\textnormal{fm}}=4$ and $\delta_{\textnormal{fm}}=8$, we can see that the classical multi-rate fixed-stress equipped with (2b) as stopping criteria needs in the last coarse mechanics step 8 iterations to converge, and between 8 and 9 iterations for the previous ones. For $\delta_{\textnormal{fm}}=2$, the classical algorithm needs 5 iterations to converge, and between 5 and 7 for the previous ones. For all the cases, the adaptive stopping criterion guarantees that the fixed-stress algorithm converged to the correct solution compared to the classical criteria (see the results with $\epsilon=10^{-3}$), while saving a substantial amount of computational effort (see the results with $\epsilon=10^{-6}$ and $10^{-4}$); see the overall performance in the three cases depicted in Figure 4(d) comparing the standard approaches with $\epsilon=10^{-6}$ to the adaptive ones.

7.1.2 Adaptive time-stepping balancing the space and time errors

In the second part of this test case, we verify the impact of the balancing criteria (4.1)-(4.2) on the fixed-stress schemes. The balancing criteria (4.1) aims adapting the times steps for the flow and mechanics subsystems in such a way that their spatial and temporal error estimators (6.26) are equilibrated through the computation. This leads practically to having $\eta^{n}_{\textnormal{sp}}\approx\eta^{n}_{\textnormal{tm}}$. That of the criteria (4.2) equilibrates the time errors from the mechanics and flow discretizations, i.e., $\eta^{n}_{\textnormal{tm,P}}\approx\eta^{n}_{\textnormal{tm,U}}$. Next, we will see that using the balancing (4.1) or (4.2) is important for the efficiency of the adaptive algorithm.

To this aim, we compare on three levels of uniform space-time mesh refinement, the standard space-time (Algo. 3.1 ) and the single- and multi-rate (Algo. 3.2) with $\delta_{\textnormal{fm}}=8,\,4$ algorithms with the adaptive fixed-stress one (Algo. 4.2). For the three refinement levels, we use the same weights $\gamma_{\textnormal{tm,J}}=0.8$ and $\Gamma_{\textnormal{tm,J}}=1.2$, J=P, U. In Figure 5 (top), the ratio of the time discretization error over the space discretization error from the flow (left) and from the mechanics (right) as a function of the total number of space–time unknowns is depicted for the aforementioned standard and adaptive algorithms. These results confirm numerically that we have distinguished the pressure and displacement errors as well as their time and space discretizations. Precisely, we can easily see that for the multi-rate schemes, the ratio $\frac{\eta_{\textnormal{tm,P}}}{\eta^{n}_{\textnormal{sp,P}}}$ (Figure 5 (top left)), remains constant when changing the ratio $\delta_{\textnormal{fm}}$, in contrast to the ratio $\frac{\eta_{\textnormal{tm,U}}}{\eta^{n}_{\textnormal{sp,U}}}$ (Figure 5 (top right)), where the ratio increases with $\delta_{\textnormal{fm}}$. The effect of the resulting ratios on the overall estimate is shown in Figure 5 (bottom left). These results make it evident that the performance of the fixed-stress algorithms is considerably improved if they are equipped with the balancing criteria (4.1) (balanced1) or (4.2) (balanced2). Particularly, the standard multi-rate algorithm reduces the computational cost of the single-rate one, but still much more expensive than the adaptive ones. In average, the adaptive one reduces the computational cost of the multi-rate one with $58\%$ while the efficiency of the algorithm in terms of precision is much more preserved. In Figure 5 (bottom right), we chose the third refinement level, then we plot the dependence of the total and fixed-stress estimators as a function of the fixed-stress iterations for the adaptive algorithms. This result confirms that the algorithm is improved if it is equipped with the balancing criteria (4.1) or (4.2). Precisely, these balancing ensure that the contribution of the fixed-stress estimator in the overall error becomes quickly negligible (see Figure 1(a) for the case without adaptivity), thus we can stop the fixed-stress iterations by setting $\gamma_{\textnormal{it}}=0.5$. Furthemore, with either of these balancing criteria we have keeping a small non-conformity in-time error which makes the application of adaptive stopping criteria more comfortable and guaranteed. Note that in the standard algorithms (as shown in the results of subsection 7.1.1), the time steps and the ratio $\delta_{\textnormal{fm}}$ are mainly based on intuition and this may induce an over-refinement in-time and may increase the nonconformity-in-time errors, affecting considerably the efficiency of the fixed-stress algorithm. In Figure 6, we plot the pressure and displacement estimators as a function of the adaptive time steps. Note that if the developed algorithm is equipped with asynchronous adaptivity in space, we can significantly reduce the total computational cost, but also the total error, as this later is dominated by the space error from the discretization of the flow subsystem.

7.2 Test problem 2: a poro-mechanical behavior of an osteonal tissue

In this test case, the poroelastic model is carried out to study the hydro-mechanical behavior of an idealized osteonal tissue. This idealized structure is a group of osteons surrounded by their cement lines and embedded in the interstitial bone matrix [34, 35]. The simplified domain presents the parts of three different osteons connected by the interstitial system (IS): a half osteon (O1) is located at the bottom of the picture and two quarters of osteons (O2) and (O3) are placed on the top-left and top-right corners, respectively (see Figure 7).

The used material properties as generated in [36] (see also [39]) are $\mathbf{K}=10^{-6}$ $(\textnormal{m}^{2})$ in the osteons and $\mathbf{K}=10^{-7}$ $(\textnormal{m}^{2})$ in the (IS)-domain. The remaining parameters are $c_{0}=0.263$ $(\textnormal{GPa})$, $\alpha=0.132$ $(\textnormal{Kg.m}^{2})$, $\mu=0.328$, and $\lambda=0.25$. The boundary conditions are $p=0$ and ${\boldsymbol{\sigma}}{\mathbf{n}}=0$ on the portion BC1 and ${\mathbf{u}}\cdot{\mathbf{n}}=0$ together with $({\boldsymbol{\sigma}}{\mathbf{n}})\cdot{\boldsymbol{\tau}}=0$, and ${\mathbf{w}}\cdot{\mathbf{n}}=0$ on BC2. The final time is $T=15~{}(\mu\textnormal{s})$.

We use Algorithm 4.2 equipped with (4.1) where we consider two computations that differ by the balancing parameters $\gamma_{\textnormal{tm,J}}$ and $\Gamma_{\textnormal{tm,J}}$. We start with an initial time step $\tau_{\textnormal{f}}^{0}=2\cdot 10^{-3}~{}(\mu\textnormal{s})$, and $\tau_{\textnormal{m}}^{0}=4\tau_{\textnormal{f}}^{0}$. The estimators are computed every 3 iterations to reduce the computational cost. Table 1 compares the number of space–time unknowns (number of asynchronous time steps, counting repetitions in the adaptive algorithm, fixed space unknowns) and performed fixed-stress iterations, and the values of the error estimators of the three computations. We observe that the gain in the number of fixed-stress iterations as well as in the number of unknowns is significant. Indeed, the two adaptive computations need approximately 30 fixed-stress to converge while the standard fixed-stress algorithm needs more than 132 iterations, thus, the total computational cost is reduced of $88.6\%$ for the first adaptive computation and of $82.5\%$ for the second one.

To clarify this gain, we can observe in Figure 8 (left) that as soon as we perform $\approx 30$ fixed-stress iterations, the fixed-stress estimator is sufficiently small to not contribute significantly on the overall error. Also as expected, the adaptive stopping criterion stops the fixed-stress algorithm when the solution is sufficiently accurate. Figure 8 (left) confirms also the role of the adaptivity in time, with which, the fixed-stress estimator becomes quickly smaller than the space and time discretization estimators, even with a large value of $\gamma_{\textnormal{it}}$, for example $\gamma_{\textnormal{it}}=0.5$. We can also observe that even with a small parameter $\gamma_{\textnormal{it}}=0.01$, the gain in the number of fixed-stress iterations is significant. In Figure 8 (right), we plot the pressure and displacement estimators as a function of time. Clearly, the displacement error dominates the pressure error along the simulation. In Figure 9, we plot the approximate solution at the final time $t=T$. Figure 10 compares the spatial discretization errors for the pressure (top left) and displacement (top right), and the fixed-stress estimator (bottom), after using our adaptive stopping criteria at the final time $t=T$. Besides detecting the dominating error at the circular boundary of the Osteons, we can see that the total error is dominated by the mechanics discretization error, and that the fixed-stress estimator is negligible.

8 Conclusion

We proposed in this paper adaptive fixed-stress iterative coupling schemes for the Biot system. Our adaptive algorithm can be used either globally-in-time or (partially) via time windowing techniques, and works as follows:

• •

At the first iteration, both time step size of flow and mechanics will be adapted in such a way that the space and time error contributions are equilibrated.

• •

We then continue iterating, where several estimators (space, time and fixed-stress) are computed, until the fixed-stress estimator becomes smaller (up to a user-chosen constant) than the other error components.

The numerical experiments demonstrated the accuracy of the estimated quantities while highlighting the applicability of the presented adaptive algorithm. Particularly, the algorithm saves important number of iterations, reduces significantly the total computational cost by adapting asynchronously the flow and mechanics time-steps and avoiding over-in-time refinement together with maintaining a small non-conformity in-time error. The algorithm may also help optimizing the tuning parameter. These benefits, together with the fact that we, a posteriori, estimate the overall error that is guaranteed and without unknown constant, leads to efficient and optimized adaptive fixed-stress coupling algorithm. Note that the present approach can be extended easily to other inexact coupling methods such as drained split, undrained split, and fixed-strain split methods. Also, the present algorithm can be applied directly, without further developments to any flux- and stress-conforming discretizations of the flow and mechanics such that cell-centered finite volume or mimitic finite difference and can easily be extended to conforming methods.

Acknowledgment

The research is supported by the Norwegian Research Council Toppforsk project 250223 (The TheMSES project: https://themses.w.uib.no).