Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem © 2020 The Authors. Published by Elsevier Ltd. Published version Kumar, Kundan; Kyas, Svetlana; Nordbotten, Jan Martin; Repin, Sergey Kumar, K., Kyas, S., Nordbotten, J. M., & Repin, S. (2021). Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem. Computers and Mathematics with Applications , 91, 122-149. https://doi.org/10.1016/j.camwa.2020.05.005 2021 Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. Computers and Mathematics with Applications xxx (xxxx) xxx Contents lists available at ScienceDirect Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa Guaranteed and computable error bounds for approximations constructed by an iterative decoupling of the Biot problem Kundan Kumara, Svetlana Kyasb, Jan Martin Nordbottena,∗, Sergey Repinc,d,e aDepartment of Mathematics, University of Bergen, Norway bGeothermal Energy and Geofluids, Institute of Geophysics, ETH Zürich, Switzerland cUniversity of Jyvaskyla, Finland dSteklov Inst. Math. of the RAS at St. Petersburg, Russia ePeter the Great Polytechnic University, St. Petersburg, Russia article info Article history: Available online xxxx Keywords: Biot problem Fixed-stress split iterative scheme a posteriori error estimates Contraction mappings abstract The paper is concerned with guaranteed a posteriori error estimates for a class of evolutionary problems related to poroelastic media governed by the quasi-static linear Biot equations. The system is decoupled by employing the fixed-stress split scheme, which leads to an iteratively solved semi-discrete system. The error bounds are derived by combining a posteriori estimates for contractive mappings with functional type error control for elliptic partial differential equations. The estimates are applicable to any approximation in the admissible functional space and are independent of the discretization method. They are fully computable, do not contain mesh-dependent constants, and provide reliable global estimates of the error measured in the energy norm. Moreover, they suggest efficient error indicators for the distribution of local errors and can be used in adaptive procedures. ©2020The Authors.Published byElsevier Ltd.This isan openaccess articleunder theCC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction The problems defined in a poroelastic medium contribute to a wide range of application areas, including simulation of oil reservoirs, prediction of environmental changes, soil subsidence and liquefaction in earthquake engineering, well stability, sand production, waste deposition, hydraulic fracturing, CO2sequestration, and understanding of the biological tissues in biomechanics. In recent years, mathematical modeling of poroelastic problems has become a highly important topic because it helps engineers to understand and predict complicated phenomena arising in such media. However, numerical schemes designed for any of the existing models provide approximations that contain errors of different nature, and these errors must be controlled. Therefore, a reliable quantitative analysis of poroelasticity problems requires efficient and computable error estimates that can be applied for various approximations and computation methods. Mathematical modeling of poroelasticity is usually based on the Biot model that consists of the quasi-static elasticity problem coupled with an equation governing slow fluid motion. Computational errors in one part of the model may seriously affect the accuracy of the other part. Therefore, getting reliable and efficient a posteriori error estimates is generally much more complicated for coupled problems than for a single equation. The Biot model is a system describing the flow and displacement in a porous medium by the momentum and mass conservation equations. Initially, it was derived at a macroscopic scale (with inertia effects negligible) in the works by ∗Corresponding author. E-mail addresses: [email protected] (K. Kumar), [email protected] (S. Kyas), [email protected] (J.M. Nordbotten), [email protected] (S. Repin). https://doi.org/10.1016/j.camwa.2020.05.005 0898-1221/©2020 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/). Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 2K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Terzaghi [1] and Biot [2]. Settlement of different types of soils was predicted in [1], which was later extended to the generalized concept of consolidation [2,3]. A comprehensive discussion of the theory of poromechanics can be found in [4]. Thus, to model the solid displacement uand the fluid pressure p, we consider the system that governs the coupling of an elastic isotropic porous medium saturated with slightly compressible viscous single-phase fluid −div(λ(divu)I+2µε(u)−αpI)=fin Q:= Ω×(0,T), ∂t(βp+αdivu)−divK∇p=gin Q,(1.1) where Qdenotes a space–time cylinder (with bounded domain Ω⊂Rd,d= {2,3}having a Lipschitz continuous boundary ∂Ωand a given time-interval (0,T), 0 <T<+∞), f∈H1(0,T;[L2(Ω)]d) and g∈L2(0,T;L2(Ω)) are the body force and the volumetric fluid source, respectively.1The first equation in (1.1) follows from the balance of linear momentum for the total Cauchy stress tensor σpor := σ(u)−αpI that accounts not only for ubut also for the pressure pscaled by the dimensionless Biot–Willis coefficient α > 0. The stress tensor for the material is governed by Hooke’s law σ(u):= 2µε(u)+λtrε(u)I=2µε(u)+λ(divu)I, where ε(u):= 1 2(∇u+(∇u)T)is the tensor of small strains. Here, λ, µ > 0 are the Lamé constants proportional to Young’s modulus Eand Poisson’s ratio νvia relations µ=E 2 (1+ν)and λ=Eν (1+ν) (1−2ν). The second equation is the fluid mass conservation (continuity) equation in Q. Here, βstands for the storage coefficient and Kis the permeability tensor assumed to be symmetric, uniformly bounded, anisotropic, heterogeneous in space, and constant in time, i.e., λK|τ|2≤K(x)τ·τ≤µK|τ|2, λK, µK>0,for all τ∈Rd.(1.2) Let Σ=∂Ω×(0,T) be a lateral surface of Q, whereas Σ0:= ∂Ω×{0}and ΣT:= ∂Ω×{T}define the bottom and the top parts of the mantel, such that ∂Q=Σ∪Σ0∪ΣT. Initial conditions are assumed to be as follows: p(x,0) =p◦∈H1(Ω) and u(x,0) =u◦∈ [H1(Ω)]don Σ0.(1.3) We introduce the following partitions of the boundary: ∂Ω=Σp D∪Σp N=Σu D∪Σu N, where Σp Dand Σu Dmust have positive measures, i.e., |Σp D|,|Σu D|>0, with the corresponding boundary conditions (BCs): p=pDon Σp D, −K∇p·n=zNon Σp N, u=uDon Σu D, σpor ·n=tNon Σu N. (1.4) For the fluid content βp+αdivu, we prescribe the following initial condition η(x,0) := βp(x,0) +αdivu(x,0) =βp◦+αdivu◦, where p◦and u◦are defined in (1.3). To simplify the exposition, we consider only homogeneous BCs, i.e., pD,zN=0 and uD,tN=0for the time being, even though all results are valid for more general assumptions. The work [5] provides the results on existence, uniqueness, and regularity theory for (1.1)–(1.4) in the Hilbert space setting, whereas [6] extends the recent results to a wider class of diffusion problems in poroelastic media with more general material deformation models. Corresponding a priori error estimates can be found in [7]. The considered system can be regarded as the singular limit of the fully dynamic Biot–Allard problem (see the details in [8]), where the acceleration of solid in the mechanics part of (1.1) is neglected. Since the Biot model is a coupled system of partial differential equations (PDEs), both iterative and monolithic approaches can be used to solve the problem (see, e.g., [9]). For the first approach, the problem can be reformulated with a contractive operator, which naturally yields iterative methods for its solution (see [8]). At each step in time, the flow problem is considered first. Next, we solve the mechanics using the pressure from the first step. The procedure is repeated until the desired convergence is reached. Different alteration of iterative cycles in flow and mechanics, i.e., single- [10] and multi-rate schemes [11,12], can be considered. The second approach is fully coupled and considers the system (1.1)–(1.4) with two unknowns simultaneously. The iterative coupling offers several advantages over the monolithic method in code design, in particular, in terms of availability of highly developed discretization methods (primal [13–15], mixed [7,16,17], Galerkin least-squares [18], finite volume (FV) [19], discontinuous Galerkin (dG) methods [20], Hybrid High-Order methods [21] isogeometric analysis [22], as well as combinations of above-mentioned ones) and algebraic solvers (e.g., general Schur complement based preconditioners [23–30] and the recently developed robust ones with respect to (w.r.t.) the model parameters 1For convenience of the reader, we collected the definitions related to the functional spaces in the Appendix. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 3 [31–35]). In the fully-coupled methods, constructing efficient preconditioning techniques for the arising algebraic systems remains a matter of active ongoing scientific research (see, e.g., [30,36–39]). The a posteriori error control for poroelastic models has been already addressed using different techniques. Application of residual-based error estimates to coupled elliptic–parabolic problems can be traced back to works [40,41]. Recently, similar error indicators were used in [42–46] for immiscible incompressible two- or multi-phase flows in porous media to address the questions of adaptive stopping criteria and mesh refinement. In [47], authors suggested an a posteriori error estimator based on a corresponding dual problem in space–time for a coupled consolidation problem that involves large deformations. In [48–51], adaptive space–time algorithms relying on the equilibrated fluxes technique were applied to the Biot’s consolidation model (formulated as a system with four unknowns). In this work, however, we turn to the functional error estimates (majorants) that are fully computable and provide guaranteed bounds of errors arising in numerical approximations. The derivation of such estimates is based on functional arguments and variational formulation of the problem in question. Therefore, the method does not use specific properties of approximations (e.g., Galerkin orthogonality) and special properties of the exact solution (e.g., high regularity). The estimates do not contain meshdependent constants and are valid for any approximation in the natural energy class. Moreover, the majorants also yield an efficient error indicator that can be used to drive mesh adaptation. Since a concise mesh adaptation algorithm is still the matter of ongoing research, we postpone including a specific example with adaptivity discussions until the next paper on this matter. Our main goal is to deduce efficient a posteriori error estimates for the approximation of the system (1.1) and demonstrate their suitability through the application to numerical problems. In [52], a posteriori error estimates of the functional type were derived for the stationary Barenblatt–Biot model of porous media. This paper deals with a more complicated Biot problem presented by an elliptic–parabolic system of partial differential equations. Our approach is based on the contraction property of the iterative method [53], which is rather general and not restricted only to the fixed-stress scheme, and functional type estimates of each equation in the Biot system (see, e.g., [54]). To the best knowledge of the authors, it is the first study targeting such a coupling between the elastic behavior of the medium and the fluid flow in the context of functional error estimates. The main result is presented in Theorem 2. Moreover, these estimates serve not only as reliable estimates of the global error measured in the energy norm, but also as efficient indicators of the local error distribution over the computational domain. The latter property makes functional majorants advantageous in automated adaptive mesh generation algorithms. The paper has the following structure: Section 2is dedicated to the generalized formulation of the Biot system and its semi-discrete counterpart derived after applying the explicit Euler scheme in time. In Section 3, we introduce an incremental approach, namely, the fixed-stress split scheme, for discretizing the considered coupled system. In particular, it clarifies the arguments for choosing the optimal parameters in the iterative scheme and proves that it is a contraction with an explicitly computable convergence rate. For the reader’s convenience, Section 4summarizes the main results of the work and the concepts that were used for their derivation, as well as highlights the most important properties of the error estimates. Sections 5and 6are dedicated to the derivation of auxiliary lemmas used in the proof of Theorem 2 (or Theorem 3) with general estimates for the approximations generated by the fully decoupled iterative approach. Finally, Section 7contains a collection of examples that illustrate the application of derived error estimates to the Biot problem. 2. Variational formulation and discretization We study approximations of the system (1.1), where ˜ V≡H1(0,T;[H1(Ω)]d) denotes the space for u(field of displacements) and ˜ W≡H1(0,T;H1(Ω)) is the space for the variable p(pressure). The generalized setting of (1.1) is read as: find a pair (u,p)∈˜ V0ט W0such that 2µ(ε(u),ε(v))Q+λ(divu,divv)Q+α(∇p,v)Q=(f,v)Q,∀v∈˜ V0,(2.1) (K∇p,∇w)Q+(∂t(βp+αdivu), w)Q=(g, w)Q,∀w∈˜ W0,(2.2) where ˜ V0:= {v∈H1(0,T;[H1(Ω)]d)|v(t)|Σu D=0a.e. t∈(0,T)}, ˜ W0:= {w∈H1(0,T;H1(Ω)) |w(t)|Σp D=0 a.e. t∈(0,T)}. The Biot system of type (2.1)–(2.2) was analyzed by several authors to establish the existence, uniqueness, and regularity of its solution. First theoretical results on the existence and uniqueness of a (weak) solution are presented in [55] for the case of β=0. Further work in this direction can be found in [5,56]. The well-posedness of the quasi-static Biot system is ensured under the above-mentioned assumptions. In fact, [57,58] established contractive results in suitable norms for the iterative coupling of (2.1)–(2.2). For an overview of the stability of existing iterative algorithms, we refer the reader to [59,60]. The system (2.1)–(2.2) can be viewed as a two-field formulation of the poroelasticity problem. In numerical analysis, there are alternative approaches that treat such a system as three- and four-field formulations. In the three-field model, an additional variable is introduced to represent the flux in the flow equations, whereas the four-field approach considers stress as yet another unknown. The three-field formulation is rather flexible since it allows different combinations of Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 4K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx discretizations. The scientific community agrees that it provides more physical approximations of the unknowns than the two-field case. Recently, the four-field formulation, where both equations were treated with mixed methods, drew much attention from the research community. The advantages of the latter representation are local conservation of mass and momentum balance and a more accurate representation of fluxes and stresses. The choice of the formulation (from the above-mentioned list) is usually motivated by the considered application as well as the limits of the computational resources. For instance, the mixed formulation of (2.2) not only provides the flux that satisfies the local mass conservation property but also generates an effective approximation of this function, which is advantageous for the functional type error control. It minimizes the majorant related to the pressure variable (see (5.2)). The same principle works for the stress field reconstruction in mechanics part. The system (2.1)–(2.2) is considered in the time-interval [0,T]divided by Nsub-intervals, such that it forms the corresponding set TN= ∪N n=1In,In=(tn−1,tn). Let un(x)∈V0and pn(x)∈W0, where V0:= {v∈V≡ [H1(Ω)]d|v⏐⏐Σu D=0}and W0:= {w∈W≡H1(Ω)|w⏐⏐Σp D=0},(2.3) respectively, are spatial parts of the solution at t=tn. Then, the semi-discrete approximation of (2.1)–(2.2) reads as (2 µε(un),ε(v))Ω+(λdivun,divv)Ω+α(∇pn,v)Ω=(fn,v)Ω,∀v∈V0, (K∇pn,∇w)Ω+1 τn(β(pn−pn−1)Ω+αdiv(un−un−1), w)Ω=(gn, w)Ω,∀w∈W0, where τn=tn−tn−1. This system generates the following problem to be solved on each step of the time-incremental method: find the pair (u,p)n∈V0×W0 (2 µε(un),ε(v))Ω+(λdivun,divv)Ω+α(∇pn,v)Ω=(fn,v)Ω,∀v∈V0,(2.4) (Kτn∇pn,∇w)Ω+β(pn, w)Ω+α(div un, w)Ω=(˜ gn, w)Ω,∀w0∈W0,(2.5) where Kτn:= τnK, the right-hand side of (2.5) is defined as ˜ gn=τngn+βpn−1+αdivun−1,(2.6) and the pair (u,p)n−1∈V0×W0is given by the previous time step. The initial values are chosen as (u,p)0=(p◦,u◦). Since, from now on, we deal only with the semi-discrete counterpart of the Biot problem, we omit the subscript Ωin the scalar product. Moreover, we always consider (2.4)–(2.5) on nth time step, which allows us to also neglect the superscript nfor the rest of the paper and consider the system (2 µε(u),ε(v)) +(λdivu,divv)+α(∇p,v)=(f,v),∀v∈V0,(2.7) (Kτ∇p,∇w)+β(p, w)+α(div u, w)=(˜ g, w),∀w∈W0.(2.8) This work aims to derive a fully guaranteed a posteriori estimates of the error between the obtained approximations (˜ u,˜ p)∈V0h×W0h, where V0hand W0hare discretization spaces of conforming approximations of functional spaces V0 and W0, respectively, and the pair of the exact solutions (u,p) of the Biot system, which is accumulated from the errors on all Ntime steps, i.e., eu:= u−˜ uand ep:= p−˜ p. On each sub-interval, these errors are measured in terms of the combined norm ⏐⏐[(eu,ep)]⏐⏐:= |||eu|||2 u+|||ep|||2 p.(2.9) In turn, each contribution is defined as follows: |||eu|||2 u:= ∥ε(eu)∥2 2µ+∥div(eu)∥2 λand |||ep|||2 p:= ∥∇ep∥2 Kτ+∥ep∥2 β,(2.10) where ∥w∥2 λ:= ∫Ωλ w2dx,∥ε(w)∥2 2µ:= ∫Ω2µε(w):ε(w) dx, and ∥w∥2 Kτ:= ∫ΩKτw·wdxare L2-norms respectively weighted with 2µ,λ, and the tensor Kτfor any scalar- and vector-valued functions wand w. The global bound of the errors euand epcontains incremental contributions from each time-interval, i.e., ∑ n=1,...,N⏐⏐[(e(n) u,e(n) p)]⏐⏐=: ⏐⏐[(eu,ep)]⏐⏐≤M(˜ u,˜ p):= ∑ n=1,...,N M(n)(˜ u,˜ p).(2.11) For the iterative approach, on each time-step In, the Biot system is decoupled into two sub-problems, where one is related to the linear elasticity, and the other — to the single-phase flow problem. Then, an iterative procedure is applied to obtain the pair (ui,pi)=(u,p)i. Next, each equation is discretized and solved, such that, instead of (u,p)i, we use the pair (u,p)i h, which contains the approximation error of the numerical method. In Section 5, we derive computable a posteriori estimates for this pair of the approximate solution. The functional M(n)is derived by combining the estimates obtained for the contractive mapping [53] and the a posteriori error majorants for elliptic problems (initially introduced [61,62]). The validity of such estimates is based on the contraction property of a specifically constructed linear combination of displacement and pressure α γdivui−L γpi,L, γ > 0, (the so-called volumetric mean stress). The selection of parameters Land γis justified and explained in Section 3. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 5 Remark 1. For the alternative monolithic approach, one solves (2.7)–(2.8) for pressure and displacement simultaneously, reconstructing the pair of approximations (˜ u,˜ p)=(uh,ph)=(u,p)h. For this case, we can derive the corresponding computable bound of the error between (˜ u,˜ p) and the exact solution. Such a functional error bound is a combination of a posteriori error estimates for each of the unknowns in (2.7) and (2.8) (see, e.g., [54,63] and references therein). Remark 2. We note that due to the Korn and Friedrichs inequalities, both ∥eu∥2and |||eu|||2 uare estimated by ∥ε(eu)∥2. Moreover, the physical bound on the Lamé parameters is given as dλ+2µ > 0, in the most general case, allowing for the first parameter λto be slightly negative for so-called auxetic materials. In this case, we use the fact that |||eu|||2 u:= ∥ε(eu)∥2 2µ+∥div(eu)∥2 λ (7.5) ≤(2 µ+dλ)∥ε(eu)∥2 holds, and work with the positively-weighted norm ∥ε(eu)∥2. However, as auxetic materials are rare, this paper will not consider such cases. Consequently, the proofs below are based on the non-negative Lamé parameters. 3. The fixed-stress splitting scheme The formal application of the iterative method to (2.7)–(2.8) yields the problem to be solved at the ith iteration step: (Kτ∇pi,∇w)+β(pi, w)+α(divui−1, w)=(˜ g, w),∀w∈W0,(3.1) (2 µε(ui),ε(v)) +(λdivui,divv)+α(∇pi,v)=(f,v),∀v∈V0,(3.2) where the flow equation (3.1) is solved for pi, using ui−1, and the elasticity equation (3.2) is used to reconstruct uiusing pirecovered on the previous step. To obtain the initial data for the iterative procedure, we first set the pressure equal to the hydrostatic pressure and then obtain u0solving (3.2). The iterative procedure proposed in (3.1)–(3.2) is known as the fixed strain, and is only conditionally stable. To stabilize the iterative scheme (3.1)–(3.2), we consider the ‘fixed-stress splitting approach’, whose properties were initially studied in [60] and [57]. This scheme operates with a special quantity: the volumetric mean total stress ηi=α γdivui−L γpi∈W,(3.3) where γand Lare certain positive tuning parameters. These parameters are usually kept constant on each half-time step. The optimal choice of γand Lproves that this iterative scheme is a contraction in the L2-norm ∥δηi∥2, where δηi:= ηi−ηi−1. Moreover, it reduces the number of iterations. By adding L(pi−pi−1) to the right-hand side of (3.1), we rewrite the system (3.1)–(3.2) using the definition (3.3) as (Kτ∇pi,∇w)+(β+L)(pi, w)=(˜ g−γ ηi−1, w),∀w∈W0,(3.4) (2µε(ui),ε(v))+(λdivui,divv)=(fi−α∇pi,v),∀v∈V0,(3.5) with complemented mixed BCs pi=0 on Σp Dand Kτ∇pi·n=0 on Σp Nas well as ui=0on Σu Dand σi por ·n=0on Σu N. Let δ∇pi:= ∇pi−∇pi−1,ε(δui):= ε(ui)−ε(ui−1), δηi:= ηi−ηi−1.(3.6) Theorem 1 establishes a contraction-type inequality for the norm ∥δηi∥2. Theorem 1 ([57,58]).If γ=α √λand L ≥α2 2λ, then the scheme (3.4)–(3.5) is a contraction that satisfies the estimate ∥ε(δui)∥2 2µ+q∥∇δpi∥2 Kτ+∥δηi∥2≤q2∥δηi−1∥2,q=L β+L,(3.7) where δ∇pi,ε(δui),δηiare defined in (3.6). Remark 3. The estimates in Theorem 1, satisfying the contraction estimate (3.7), also hold for the Galerkin approximations {δηh}i∈Wh, where Whis a discretization space of W. Moreover, in the Appendix, we show that a similar contraction theorem holds for the sequence {δ(η−ηh)i} ∈ Wh, where ηi∈Wis the generated by the fixed-stress split iterative scheme defined in (3.4)–(3.5) and ηi h∈Whis discretization of the latter sequence. Generally, it is important to note that all theorems and lemmas below are formulated for a pair (u,p)i∈V0×W0that forms a contraction w.r.t. to (u,p)i−1∈ V0×W0and its discrete approximation (u,p)i h∈V0h×W0hthat forms a contraction relative to (u,p)i−1 h∈V0h×W0h. Remark 4. There exist alternative ways to choose the tuning parameter L. In particular, the physically motivated choice Lcl =α2 λ+2µ/dis considered in [60], whereas [58] suggests Lopt =α2 2 (λ+2µ/d). The recent study [64] suggests the numerical evidence on the iteration counts w.r.t. the full range of the Lamé parameters for heterogeneous media. Numerical investigation of the optimality of these parameters and their comparison with physically and mathematically motivated values from the literature was done in [65]. The authors demonstrated that their optimal value is dependent not only on mechanical material parameters but also on the boundary conditions and material parameters associated with the fluid flow problem. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 6K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Remark 5. The inequality (3.7) shows that the sequence {δηi}i∈Nis generated by a contractive operator. Therefore, due to the Banach theorem, it tends to a certain fixed point. Moreover, since all terms on the left-hand side of (3.7) are positive, in practice, {δηi}i∈Nmight converge with an even better contraction rate than q=L β+L. Corollary 1. From Theorem 1, it follows that ∇δpiand ε(δui)in (3.6) are also converging sequences and satisfy ∥∇δpi∥2 Kτ≤q∥δηi−1∥2and ∥ε(δui)∥2 2µ≤q2∥δηi−1∥2, respectively. We use Corollary 1 to derive the error estimate for the term |||ep|||2 p. In particular, it yields the following result based on the estimates for the Banach contractive mappings (see [53,54]). Lemma 1 (Estimates for Contractive Mapping).Let the assumptions of Theorem 1hold. Then, we have the estimates ∥∇(p−pi)∥2 Kτ≤q (1−q)2∥δηi−1∥2,(3.8) ∥ε(u−ui)∥2µ≤q2 (1−q)2∥δηi−1∥2.(3.9) Proof. Consider ∥∇(pi+m−pi)∥Kτ≤ ∥∇(pi+m−pi+m−1)∥Kτ+. . . +∥∇(pi+1−pi)∥Kτ ≤q(∥ηi+m−1−ηi+m−2∥+. . . +∥ηi−ηi−1∥) ≤q(qm+. . . +1)∥ηi−ηi−1∥. By taking the limit m→ ∞ and noting that in this case (qm+qm−1+···+1) →1 1−q, we arrive at (3.8). The inequality (3.9) is proved using similar arguments. □ If in Lemma 1 we consider the iterations iand i−m,i>mas two subsequent iterations, a more general version of the estimates (3.8) and (3.9) can be formulated. Lemma 2 (General Estimates for Contractive Mapping).Let the assumptions of Theorem 1hold. Then, we have the estimates ∥∇(p−pi)∥2 Kτ≤min 1≤m≤i{qm (1−qm)2∥ηi−ηi−m∥2},(3.10) ∥ε(u−ui)∥2µ≤min 1≤m≤i{q2m (1−qm)2∥ηi−ηi−m∥2}.(3.11) Proof. Proof follows along the lines of the proof of Lemma 1 with m=1,...,i.□ Remark 6. The estimates in Lemma 2 improve the value of qm (1−qm)2and q2m (1−qm)2if qis close to 1. At the same time, it might look counterintuitive, but the choice m=iis not always optimal. As the quotient with qdecreases, the term ∥ηi−η0∥2 grows. Therefore, we should choose mcarefully. Computing the majorant on each step of our iterative algorithm might be computationally expensive, therefore, in certain cases, mmust be chosen a priori. Alternatively, with several extra iterations after reaching the desired convergence in pand u, both q′ (1−q′)2with q′ = qmand ∥ηi−ηi−m∥2will decrease, impacting total values of the majorant. Remark 7. Theorem 1 also yields the so-called a priori contractive estimates, i.e., ∥∇(p−pi)∥2 Kτ≤q2i−1 (1−q)2∥η1−η0∥2,(3.12) ∥ε(u−ui)∥2µ≤q2i (1−q)2∥η1−η0∥2,(3.13) which can be used as an alternative upper bound. 4. Main results Theorem 2 presents the main result of this work, that is an upper bound of the error ⏐⏐[(eu,ep)]⏐⏐(cf. (2.11)). Theorem 2 (On Functional Error Estimates).For any pi h∈W0and ui h∈V0, we have the estimates ∑ n=1,...,N⏐⏐[(e(n) u,e(n) p)]⏐⏐=: ⏐⏐[ep,eu]⏐⏐≤M:= ∑ n=1,...,N M(n) p+M(n) u, Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 7 where M(n) pand M(n) ucorrespond to the error introduced by approximation schemes in each variable p and u, respectively, and measured by the norms |||e(n) u|||2 uand |||e(n) p|||2 p. Omitting the suffix (n)for readability and generality, we can formulate both functionals as Mp:= 2(Mh p(pi h,zi h)+min {Mi p,Mi,m p,M : i p}), Mu:= 2(Mh u((u,p)i h,(τ,z)i h)+min {Mi u,Mi,m u,M : i u}), where Mh p,Mi p,Mi,m p, and M : i pare defined in Lemmas 3,6,Corollary 4, and Lemma 7, respectively, as well as Mh u,Mi u,Mi,m uand M : i uare presented by Lemmas 5,8,Corollary 5, and Lemma 9, respectively. The functionals Mh pand Mh uare defined in Section 5by means of functional arguments. They provide the upper bounds of errors introduced in (3.4)–(3.5) for the ith iteration, when the system is solved numerically. To be precise, these error functionals provide a bound between the exact solution (u,p)i=(ui,pi)of(3.4)–(3.5) and its approximation (u,p)i h=(ui h,pi h). Mh pand Mh uare reliable and dependent only on explicitly computable constants, approximations, and auxiliary functions. Estimates Mi p, Mi,m p, and M : i p, as well as Mh u, Mi u, and M : i uare derived in Section 6by means of the contraction mappings estimates [53]. Bounds Mi pand Mi ufollow naturally from Lemma 1, whereas Mi,m pand Mi,m uare derived from a more general Lemma 2. Derivation of the functionals M : i uand M : i uis based on yet another basic property of the contractive operators (3.12)–(3.13) that is highlighted in Remark 7. Generally, Theorem 2 provides a mathematical tool for reliable error control of approximate solutions of the Biot problem in the poroelastic medium. The functional M provides a guaranteed bound of the error in these approximations, which is confirmed by numerical examples in Section 7. The set of computational tests is designed to provide an overview of several important properties of functional error estimates, as well as confirm the universality w.r.t. some of the parameters coming either from the mathematical model (e.g., permeability, Lame coefficients, etc.) or dictated by the iterative scheme (tuning parameters of the fixed-stress split method). 5. Estimates of errors generated by discretization Before deriving the estimates of approximation errors that appear in the contractive iterative scheme, we need to study the discretization errors encompassed in (3.4)–(3.5) for the ith iteration. Henceforth, the pair (u,p)i=(ui,pi) is considered as the exact solution of (3.4)–(3.5), whereas (u,p)i h=(ui h,pi h) denotes its approximation computed by a certain discretization method. We aim to derive computable and reliable estimates of the error measured in the terms |||ei p|||2 pand |||ei u|||2 u. Majorant of the error in the pressure term. For the first equation (3.4),Lemma 3 presents a computable upper bound of the difference ei p:= pi−pi h between the exact solution pi∈W0and its approximation pi h∈W0, measured in terms of the energy norm |||ei p|||2 p. Lemma 3. For any pi h∈W0, any auxiliary vector-valued function zi h∈HΣp N(Ω,div) := {zi h∈ [L(Ω)]d|divzi h∈L2(Ω),zi h·n∈L2(Σp N)},(5.1) and any parameter ζ≥0, we have the following estimate ∥∇ei p∥2 Kτ+∥ei p∥2 β=: |||ei p|||2 p≤Mh p(pi h,zi h;ζ), where Mh p(pi h,zi h;ζ):= (1 +ζ)∥rd(pi h,zi h)∥2 K−1 τ+(1 +1 ζ)Cp Ω(∥req(pi h,zi h)∥2 Ω+Mh q((u,p)1 h)+∥zi h·n∥2 Σp N).(5.2) Here, rd(pi h,zi h):= zi h−Kτ∇pi h,req(pi h,zi h):=˜ g−γ ηi−1 h−(β+L)pi h+divzi h, where ˜ g is defined in (2.6), and Mh q:= (Cq(α γMh,1/2 p,L2+L γMh,1/2 u,div )+(Cq+1) ∥η0−η0 h∥)2 ,Cq:= i−1 ∑ k=1 qk+1, Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 8K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx where Mh p,L2(p1 h,z1 h)and Mh u,div((u,p)1 h,τ1 h,z1 h), defined in Corollaries 2and 3, are dependent on the explicitly given η0. The constant (Cp Ω)2:= 1 (β+L)(1+(Ctr Σp N)2)(5.3) is defined via the constant in the trace-type inequality ∥w∥Σp N≤Ctr Σp N∥w∥Ω,∀w∈W0,(5.4) and the positive parameters of the Biot model, βand L. Proof. The majorant Mp(pi h,zi h;ζ) follows from [66, Section 2] and [54, Section 4.2–4.3], i.e., we consider (3.4) with a bilinear form (Kτ∇pi h,∇w)+(β+L)(pih, w) subtracted from its left- and right-hand sides (Kτ∇ei p,∇w)+(β+L) (ei p, w)=(˜ g−γ ηi−1−(β+L)pih, w)−(Kτ∇pi h,∇w). Next, we set w=ei pand introduce an auxiliary function zi h∈HΣp N(Ω,div) (cf. (5.1)) satisfying the identity (divzi h,w)Ω+ (zi h,∇w)Ω=(zi h·n,w)Σp N, such that ∥∇ei p∥2 Kτ+∥ei p∥2 β+L=(zi h−Kτ∇pi h,∇ei p)+(˜ g−γ ηi−1−(β+L)pih+divzi h,ei p)−(zi h·n,ei p)Σp N =(rd(pi h,zi h),∇ei p)+( ˜ req(pi h,zi h),ei p)−(zi h·n,ei p)Σp N,(5.5) where˜ req(pi h,zi h):=˜ g−γ ηi−1−(β+L)pi h+divzi h. Using the Hölder and Young inequalities, we can estimate the first term on the right-hand side of (5.5) as (rd(pi h,zi h),∇ei p)≤1 2(1 +ζ)∥rd(pi h,zi h)∥2 K−1 τ+1 2(1+ζ)∥∇ei p∥2 Kτ.(5.6) The second term on the right-hand side of (5.5) is bounded analogously, i.e., ( ˜ req(pi h,zi h),ei p)−(zi h·n,ei p)Σp N≤1 2(1 +1 ζ) (Cp Ω)2(∥˜ req(pi h,zi h)∥2+∥zi h·n∥2 Σp N )+1 2 ζ 1+ζ∥ei p∥2 β+L,(5.7) where Cp Ω(cf. (5.3)) is a constant in the inequality ∥w∥2+∥w∥2 Σp N≤(Cp Ω)2∥w∥2 β+L,∀w∈W0, defined in (5.4). By summing up the results of (5.6) and (5.7), we obtain ∥∇ei p∥2 Kτ+∥ei p∥2 β≤ ∥∇ei p∥2 Kτ+∥ei p∥2 β+L ≤(1 +ζ)∥rd(pi h,zi h)∥2 K−1 τ+(1 +1 ζ) (Cp Ω)2(∥˜ req(pi h,zi h)∥2+∥zi h·n∥2 Σp N ).(5.8) At this point, the term ∥˜ req(pi h,zi h)∥2is not fully computable in the usual sense of functional majorants since it is defined using ηi−1. However, it can be estimated by ∥˜ req(pi h,zi h)∥2≤ ∥˜ g−γ ηi−1−(β+L)pi h+divzi h∥2 ≤2∥˜ g−γ ηi−1 h−(β+L)pi h+divzi h∥2+2γ2∥ηi−1−ηi−1 h∥2 ≤2∥req(pi h,zi h)∥2+2γ2∥ηi−1−ηi−1 h∥2.(5.9) Consider the norm ∥ηi−1−ηi−1 h∥(without the squares) and apply an approach similar to the one that was used to prove Lemma 1 ∥ηi−ηi h∥ ≤ i−1 ∑ k=1∥δηk−δηk h∥+∥η0−η0 h∥Theorem 4 ≤(qi−1+. . . +1) ∥δη1−δη1 h∥+∥η0−η0 h∥ ≤(i−1 ∑ k=1 qk+1)∥δη1−δη1 h∥+∥η0−η0 h∥ ≤(i−1 ∑ k=1 qk+1)(∥η1−η1 h∥+∥η0−η0 h∥)+∥η0−η0 h∥ ≤Cq∥η1−η1 h∥+(Cq+1) ∥η0−η0 h∥. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 15 Proof. To decompose the error |||ep|||2 pinto two parts, we apply the triangle inequality |||ep|||2 p= |||p−pi h|||2 p≤2(|||p−pi|||2 p+|||pi−pi h|||2 p).(6.18) The first term on the right-hand side of (6.18) is bounded by (6.1) from Lemma 6, whereas the second term is controlled by (5.2) from Lemma 3. Analogously, using the triangle rule, we obtain |||eu|||2 u= |||u−ui h|||2 u≤2(|||u−ui|||2 u+|||ui−ui h|||2 u).(6.19) The first term on the right-hand side of (6.19) is controlled by (6.12),(6.16), or (6.17), whereas the estimate of the second term follows from (5.17).□ 7. Numerical examples Numerical properties of the estimates above are explained in the following series of examples. We consider three tests with manufactured solutions for pressure and displacement. First two tests correspond to different values of the contraction parameter qin the fixed-stress scheme and study the efficiency of error estimators. For both of them, we consider the material properties (permeability, the Lame coefficients, etc.) that are more academic as well as more realistic. The third test (with non-polynomial exact solutions) is considered in order to exclude the effect of the super-convergence while testing the numerical scheme. We solve the Biot model at each moment in time, using the fixed-stress scheme by choosing different mesh sizes and time steps. Moreover, we vary finite element pairs used for solving the variational problem. All these alterations, including material properties, discretization parameters, and the finite element pairs, allow us to study the performance of the estimators and show their robustness. Besides, the first and third examples include the comparison of the CPU cost needed to solve (3.4)–(3.5) for pressure and displacement and of the time effort spend on the reliable error control for each variable. In each test below, we study the convergence of the numerical scheme applied to the solver (3.4)–(3.5) for chosen discretization time step and spatial mesh size. The errors are computed by employing different norms, the total error norm, including the pair of unknowns, as well as the individual ones in both L2and the energy norms. We fix the number of fixed-stress iterations to study the convergence. The majorants Mh pand Mh uare either directly computed or minimized w.r.t. to the auxiliary functions. To characterize the efficiency of all above-mentioned error estimators, we use a so-called efficiency index defined as Ieff := M/|||e|||2, where M is a chosen majorant and |||e|||2is the error measured in the corresponding norm. We study individual contributions of different terms to the majorants, namely of the dual term and the reliability/equilibration term. We demonstrate that the first one provides quantitatively efficient error indicators, whereas the small contribution of the second one ensures the reliability of total error bounds. We also highlight the contributor of iterative majorants Mi pand Mi,m pinto the total error bounds for various ranges of contraction parameter q. In particular, the second example considers a set of parameters resulting in qclose to 1 and highlights the quantitative improvement of the total error bound when Mi,m pis used instead of Mi p. For the readers’ convenience, the time steps below are stated in seconds [s], and the spatial discretization steps are measured in meters [m]. We also assume that Ω:= (0,1)2∈R2,T=10 throughout all the examples. Example 1. Verification of majorants’ properties on simple manufactured solution (w.r.t. two different sets of parameters). First, we consider an example with polynomial manufactured solution, in which we clarify the process of both discretization and iterative error estimates calculation as well as highlight the most important properties of the introduced majorants. Therefore, we start with relatively simple material parameters that will be changed at the end of the example to test more realistic scenarios. Let the exact solution of (1.1) be defined as u(x,y,t):= t x (1 −x)y(1 −y)[1,1]Tand p(x,y,t):= t x (1 −x)y(1 −y). The Lame parameters are µ=λvol =1, which leads to λplane =2λvol µ λvol+2µ=2 3. We set α=β=1, CF=1 √2π, and K[mD/cP] is a unit tensor. From the parameters above, it follows that L=α2 2 (λ+2µ/d)=0.3 and q=L β+L=3 13 ≈0.23077. Taking into account that the error estimates depend on the term q2 (1−q2), which goes to infinity as qgoes to 1, the efficiency of the resulting Muand Mpdrastically depends on the value of q. In this particular example, the ratio q2 (1−q2)=0.05625 is relatively small, which prevents Mi uand Mi poverestimating the errors. We start with discretization by 10 time steps (where the total number of steps in time is denoted by N) with corresponding size of the step τ= 1.0 and spatial mesh-size h=1 64 using standard P1(polynomial/Lagrangian first-order) finite elements (see (7.7)). Let Idenote the number of iterations to solve (3.4)–(3.5) on each time-step; we consider I=5 for this discretization. Table 1 illustrates the convergence of the errors in uand pw.r.t. the iteration steps i=1,...,I for the time-interval [t9,t10]=[9.0,10.0]. We note here, that the study of the dependence of the iteration number on discretization parameters is beyond the scope of this paper. Therefore, in all discussed numerical tests, Iis constant. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 16 K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Table 1 Example 1 (Academic parameters). Errors and majorants at the time step [t9,t10],τ=1.0, for h=1 64 , and I=5. The values are measured w.r.t. the increment in ∥p∥2 pand ∥u∥2 u. i=1,...,I∥ep∥2Mh p∥ep∥2 βMh p,L2∥eu∥2Mh u∥diveu∥2 λMh u,div 2 1.87e−04 1.87e−04 4.77e−09 8.90e−06 1.86e−04 4.95e−04 4.45e−05 1.06e−04 4 1.87e−04 1.87e−04 2.04e−09 8.90e−06 1.86e−04 4.95e−04 4.45e−05 1.06e−04 5 1.87e−04 1.87e−04 2.04e−09 8.90e−06 1.86e−04 4.95e−04 4.45e−05 1.06e−04 Table 2 Example 1 (Academic parameters). Decrease of the values of ζand Mh p, as well as the terms m2 dand m2 eq that the majorant contains, w.r.t. the number of optimization cycles. The values of Mh pand Mh uare obtained by minimizing each functional w.r.t. to the auxiliary functions. For instance, the functional estimate to control the error in the variable p Mh p(pi h,zi h;ζ):= (1 +ζ)∥rd(pi h,zi h)∥2 K−1 τ+(1 +1 ζ)Cp Ω(2∥req(pi h,zi h)∥2 Ω+2γ2q2(i−1) ∥η0−η0 h∥2+∥zi h·n∥2 Σp N)(7.1) is minimized w.r.t. the function zi hand parameter ζ. The results of this optimization procedure are listed in Table 2. We see that only two iterations are enough to achieve a good efficiency Ieff := Mh p/|||ep|||2=1.0023. Table 2 illustrates the decrease of the dual term m2 d:= ∥rd(pi h,zi h)∥2 K−1 τ= ∥zi h−Kτ∇pi h∥2 K−1 τ and the reliability/equilibration term of the majorant m2 eq := ∥req(pi h,zi h)∥2 Ω= ∥˜ g−γ ηi−1 h−(β+L)pi h+divzi h∥2 Ω. In order to achieve the desired efficiency of the error estimate, the reliability term must be several orders of magnitude smaller than the dual one, which holds in this case. Instead of minimizing the functional Mh u, we apply a post-processing of the function ¯ ui h:= P(ui h), and then substitute it to approximate the auxiliary function τi h=Lε(¯ ui h). This yields a sufficiently efficiency index in the range 2 −6, while saving computational effort otherwise associated with calculating the minimizer of Mh u. In this example, we observe that the terms m2 d,Kand m2 d,λ,µ become rather close to the true errors |||ep|||2and |||eu|||2, respectively, and can be used as error indicators for the local error distribution over the computational domain. To confirm that, Fig. 1 presents the distribution of errors and indicators (generated by majorants) w.r.t. numbered finite element cells. On the right side of Fig. 1, we depict |||ep|||2and m2 d,K:= ∥rd,K(pi h,zi h)∥2and on the left side |||eu|||2and the indicator m2 d,µ,λ := ∥rd,µ,λ(ui h,τi h)∥2w.r.t. numbered finite element cells. Green marker depicts the error and red color represents the majorant values. We see that the indicator produces a quantitatively efficient error. First, let us assume that the error bounds Mh pand Mh p,L2, as well as Mh uand Mh u,div are reconstructed only at the 4th and 5th steps to calculate Mi uand Mi pfor i=5 (using (6.1) and (6.12)). Let |||e(n) p|||2and |||e(n) u|||2,n=1,...,Ncorrespond to the error increments at the nth time-step. The contribution of the majorants for discretization errors and the iterative majorant at the Nth time step is of the same magnitude, i.e., Mh,(N) p=1.87e−04,Mi,(N) p=3.05e−04,Mh,(N) u=4.95e−04,Mi,(N) u=3.20e−05. Then, the Nth increment of both total errors and the corresponding majorants for pressure and displacement, respectively, is as follows: |||e(N) p|||2=4.86e−05,M(N) p=2.56e−04,Ieff(M(N) p)=2.29 and |||e(N) u|||2=4.85e−05,M(N) u=2.74e−04,Ieff(M(N) u)=2.38. This makes the general error, majorant, and the corresponding efficiency index contributions to be the following [(e(N) u,e(N) p)]2=4.85e−04,M(N)=9.99e−05,Ieff(M(N))=2.36. If instead of the 4th and 5th iterations we use 2nd and 5th ones to compute Mi,m,(N) uand Mi,m,(N) p(using m=3 in (3.10) and (3.11)). Then contribution of iterative majorants is Mi,m,(N) p=5.03e−06 and Mi,m,(N) u=2.81e−08, and the values of the relative errors and majorants accumulated over all time steps are |||ep|||2=1.87e−04,|||eu|||2=1.86e−04,Mp=3.84e−04,Mu=9.91e−04, Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 17 Fig. 1. Example 1(Academic parameters). (a), (c) Local distribution of the error |||ep|||2in red and the error indicator m2 d,Kgenerated by the majorant Mh pin green, and (b), (d) local distribution of the error |||eu|||2in red and the error indicator m2 d,µ,λ generated by the majorant Mh uin green. The distributions are depicted w.r.t. numbered finite element cells of the uniform meshes with (a)–(b) h=1 /32 and (c)–(d) h=1 /64. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) respectively. The corresponding efficiency indices over the entire time interval are summarized in Table 3(a). The results are presented w.r.t. meshes with two different mesh-sizes, hand τ. As the caption suggests, all these values are obtained for auxiliary functions reconstructed by the following finite elements zi h∈RT1and τi h∈ [P2]2×2, which is an order higher than the usual choice of finite element approximation spaces for the fluxes and stresses in mixed formulations. To stress the importance of this to obtain more accurate error bounds, Table 4 lists the results obtained for cases when different finite element pairs for auxiliary functions are used, i.e., (a) zi h∈RT0and τi h∈ [P2]2×2or (b) zi h∈RT1and τi h∈ [P1]2×2. Table 5 illustrates the CPU time needed for various stages of the problem solution, i.e., first computation of both pressure and displacement, as well as a posteriori error control for each of the unknowns. For each operation, we also present the time needed as percentage of the total time spent on each time step. First, we consider Table 5(a). As expected, computing pis a much cheaper problem to solve (75–85 s on average) in comparison with the computation of u(260– 270 s). Since a direct calculation of Mh pdoes not involve the majorant minimization, the cost of error control for pis about 13–14 s. Reconstruction Mh uis more time-consuming than the one for the pressure, i.e., average time is 100–120 s. If we allow one iteration to minimize the majorant Mh p, the cost for the error control of pincreases (see Table 5(b)). The third column of the table indicates that one such iteration costs additional 200 s, and the error control for p-variable accounts for about 30%–35% of the total cost. Computational time for the case with two iterations minimizing Mh pis illustrated in Table 5(c). Here, additional 200 s are added to the third column, making it the most costly part of the entire time step. Finally, if we allow one iteration of the Mh uminimization, the cost of error control increases even more. Eventually, as we improve the efficiency index of the error bounds, the computational costs to reconstruct functionals Mh pand Mh udominates over the time needed for the calculation of Mi pand Mi u, accounting for 95% of the total error control. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 18 K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Table 3 Example 1 (Academic parameters). Convergence of errors and majorants w.r.t. the choice of spatial mesh sizes hand time steps τ(measured relatively w.r.t. ∥p∥2 pand ∥u∥2 u). For all cases, the auxiliary functions are reconstructed by the following finite elements zi h∈RT1and τi h∈ [P2]2×2. Table 4 Example 1 (Academic parameters). Convergence of errors and majorants w.r.t. the choice of approximation spaces for zi hand τi h(measured relatively w.r.t. ∥p∥2 pand ∥u∥2 u). For all cases, simulation is performed for the discretization with τ=1.0. Now, we consider N=103and h=1 64 , in particular, the last time-step [t999,t1000]of a size τ=0.01. After ten iterative steps (I=10), the error bounds are Mh,(N) p=6.66e−05,Mi,(N) p=3.04e−02,Mh,(N) u=4.90e−04,Mi,(N) u=3.15e−05. We emphasize that contributions of majorant for discretization error and the iterative majorant for the pressure variable substantially differ in magnitude. These are the results obtained considering two last subsequent iteration steps with q=0.2307 that correspond to 3q 1−q2((CF Σp D)2β/(λKτ)+1)=7.09. Such a difference in magnitudes of majorants for p results in the efficiency index Ieff(M(N) p)=43.10. However, we can exploit the flexibility of Lemma 8 with M=7 to obtain a better contraction parameter ˜ q:= q7=3.48e−05, i.e., Mi,m,(N) p=1.41e−06 and Mi,m,(N) u=2.21e−13. These values improve the efficiency index of, for example, M(N) pby approximately 21.12 times, and yield |||e(N) p|||2=9.83e−08,M(N) p=4.07e−07,Ieff(M(N) p)=2.04 and |||e(N) u|||2=5.59e−09,M(N) u=2.93e−08,Ieff(M(N) u)=2.29. The accumulated values over the whole time interval are as follows: |||ep|||2=3.28e−05,Mp=1.36e−04,and Ieff(Mp)=2.04, Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 19 Table 5 Example 1 (Academic parameters). Comparison of the CPU time required to perform different stages of numerical calculations w.r.t. the selected time steps for the discretization parameters N=10, τ=1.0, h=1 64 . Table 6 Example 1 (Realistic parameters). Convergence of errors and majorants w.r.t. the choice of spatial mesh sizes (measured relatively w.r.t. ∥pi h∥2 pand ∥ui h∥2 u) for N=100 and τ=0.1. |||eu|||2=1.86e−06,Mu=9.80e−06,and Ieff(Mu)=2.29. This yields the total values ⏐⏐[(eu,ep)]⏐⏐=2.40e−06 and M =1.2029e−05 with the efficiency index Ieff =2.24. The latter values amount others are included in Table 3(c). We see that even with a decreasing τ, which scales the permeability tensor K, the efficiency of total error estimates stays rather robust. We note here that each increment of M(N) pand M(N) ucan be computed using M : i pand M : i u, which does not require the majorant of discretization errors Mh,(N) pand Mh,(N) uon iteration steps i−1 or i−mand, as a consequence, their minimization w.r.t. the auxiliary functions zi hand τi h, respectively. Moreover, the values M : i pand M : i uare orders of magnitude smaller than Mi,m,(N) pand Mi,m,(N) u, i.e., M : i p=7.07e−10 and M : i u=7.33e−13. Such values automatically minimize the contribution of iterative estimates in both M(N) pand M(N) u, as well as further improve their efficiency. Let us assume now more realistic parameters in the example above (similar to those taken from [64] and [58]). Let the exact pressure be scaled as follows: p(x,y,t)=108x(1 −x)y(1 −y)t. The permeability tensor divided by the fluid viscosity is taken as K=100I[mD/cP], the fluid compressibility fixed to 4.7·10−7[psi−1], and initial porosity φ0is assumed to be 0.2. The Biot and the bulk modulus are M=1.65 ·1010 [Pa] and E=0.59 ·109[Pa], respectively. That results in β=1 M+cfφ0=9.40 ·10−8, which, in turn, yields a considerably small L=α2 2 (λ+2µ/d)=1.35 ·10−9. Such a tuning parameter generates instantaneous convergence of an iterative scheme with the contractive parameter q=L β+L=6.73 ·10−12. The resulting errors and corresponding estimates are summarized in Table 6 (for τ=0.1 and difference between spatial mesh-sizes h). For such q, even one iteration is enough for convergence. However, one can consider two/five iterations to improve the sharpness of the majorant. The efficiency indices obtained here confirm the quantitative properties of total majorants also for parameters close to those used in engineering applications. Example 2. Dependence of the total error bound on the iterative majorants contributions (w.r.t. two different sets of parameters). In the next example, we consider another polynomial exact solution. However, the first set of material Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 20 K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Table 7 Example 2 (Academic parameters). Errors and majorants w.r.t. the iteration steps for N=10, τ=1.0, h=1 64 , and I=12 (both values are measured relative to the increment in ∥p∥2 pand ∥u∥2 uat the Nth time step). Table 8 Example 2 (Academic parameters). Convergence of errors and majorants w.r.t. the choice of spatial mesh sizes and time steps (measured relatively w.r.t. ∥pi h∥2 p,∥ui h∥2 u, and [(eu,ep)]2). parameters (which can be regarded as more academic ones) are chosen in such a way that parameter qtakes a value close 1. Such a contraction can cause the deteriorated quality of iterative error estimates, which compromises the efficiency index of total majorants. Following similar arguments made in Example 1, we highlight the improvement in the efficiency index of total error bound M when the functional Mi,mis used instead of Mi. In particular, it is done for different time steps and spatial mesh-sizes (see Table 7). Moreover, in the second part of this example, we choose more realistic material parameters, in particular, anisotropic permeability tensor of relatively small magnitude, which might cause a quality decrease in the discretization majorants (see, e.g., [69]). Having said that, the more realistic set of parameters often lead to rather small contraction values, making iterative component of the error bounds rather neglectable. Regardless of chosen parameters, we present the efficiency indices that confirm the robustness of total error bound. Let the exact solution of (1.1) be defined as u(x,y,t):= [t(x2+y2) t(x+y)]and p(x,y,t):= t x (1 −x)y(1 −y). We fix the Poisson ratio to be ν=0.2 and Young modulus E=0.594 [Pa], which yields the Lame parameters µ=0.25 and λ=0.12. We set α=1 and β=1 M+cfφ0=0.11, where M=1.65 ·10−10 [Pa], Ku=K+α2 c0=10.28, cf=1 c0λK K+4/3µ Ku+4/3µ=0.58, CF=1 √2π, and K=I[mD/cP], where Iis a unit tensor. From the parameters above, it follows that L=α2 2 (λ+2µ/d)=1.34 and q=L β+L=0.92. With such q, the ratio q2 (1−q)2is 133.33, which might influence the quantitative performance of the majorant. We consider 10 time steps of the length τ= 1.0 and a spatial mesh-size h=1 64 using standard P1finite elements (see (7.7)). The number of iterations to solve the problem on each time step is set to 12. For the time interval [t9,t10], the convergence of errors in uand pis presented in Table 7. From one side, we can consider the iterations Iand I−1 as subsequent ones with a contraction parameter q=0.92. From the other side, by using Lemma 2 instead, let m=6 so that the 6th and 12th iterations are treated as two consecutive steps with q′=q6=0.60. Then, the constants dependent on the ratio q2 (1−q2)=2.30 attain more acceptable values. Table 8 illustrates the improved efficiency indices of error majorant as the method explained-above is employed. Local distribution of the error in pand uon each cell of the finite-element discretization is presented for mesh sizes h=1 /4and h=1 /8in Fig. 2. One can see the resemblance in the local error distribution for psince the exact solutions for both examples are the same. Local values of the error and indicator in uhave a more uniform distribution. Let us consider more realistic parameters similar to those chosen in [70]. Again, let the exact pressure be p(x,y,t):= 108x(1−x)y(1−y)t. Mechanical parameters such as the Biot and the bulk modulus are chosen as follows: M=1.45·104 Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 21 Fig. 2. Example 2(Academic parameters). Errors |||ep|||2and |||eu|||2(in red) and error indicators m2 d,Kand m2 d,µ,λ (in green) generated by the majorants Mh pand Mh u, respectively, distributions w.r.t. numbered finite element cells. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) [Pa] and E=7·107[Pa], respectively. The permeability tensor is K=diag{50,200} · 10−10 [mD/cP], whereas the fluid viscosity is µf=10−3. Thus, we obtain β≈1 M+cfφ0=6.89 ·10−5, which, in turn, yields a considerably small L=α2 2 (λ+2µ/d)=1.14 ·10−8. Such a parameter generates contractive parameter q=L β+L=1.66 ·10−4. The latter means that the ratio q2 (1−q)2is of order 10−8, which does not degenerate the efficiency of majorant. We summarize the convergence results in Table 9. It again confirms that even in case of more realistic parameters common for engineering applications, the value of estimates remains rather efficient. Example 3. Verification of the majorants properties w.r.t. non-polynomial manufactured solutions. To make sure that we exclude the super-convergence in testing the scheme presented above, we consider the non-polynomial manufactured solution. The chosen pand uin (1.1) are u(x,y,t):= [xsin πx(1 −y) sin πy t x(1 −x)y(1 −y)(sin t+1)]and p(x,y,t):= sin πxsin πy(t2+t+1). The Lame parameters are similar to those considered in Example 1, i.e., µ=λvol =1, λplane =2 3. Besides, α=β=1, CF=1 √2π, and K=I[mD/cP]. The parameters above yield L=0.3 and q=3 13 . Table 10 presents convergence results corresponding to meshes with different mesh-sizes hand time steps τ. The number of iterations is fixed to be I=5 throughout all tests. Learning from the experience of first two examples, we apply Lemma 2 with m=3, which yields a smaller contractive parameter ˜ q:= q3=0.0123. The table is divided into Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 22 K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Table 9 Example 2 (Realistic parameters). Convergence of errors and majorants w.r.t. the choice of spatial mesh sizes and time steps (measured relatively w.r.t. ∥pi h∥2 p,∥ui h∥2 u, and [(eu,ep)]2). Table 10 Example 3. Convergence of errors and majorants w.r.t. the choice of spatial mesh sizes hand time steps τ(measured relatively w.r.t. ∥p∥2 p,∥u∥2 u, and [(eu,ep)]2). two parts: Table 10(a) presents the results obtained for the auxiliary functions zi h∈RT1and τi h∈ [P2]2×2, whereas Table 10(b) uses the functions zi h∈RT0and τi h∈ [P2]1×1. As expected, the efficiency indices in the first part of the table are slightly better than in the second one, since the auxiliary functions that minimize the majorants of discretization errors are reconstructed more accurately. Besides the robustness of error estimates w.r.t. different discretization parameters, we address the computational time for various scenarios. We consider discretization time steps τ=1 (see Table 11(a)) and τ=0.1 (Table 11(b)). In the first part of the table, corresponding to a smaller step, we first illustrate the CPU costs for a more computationally-heavy error control, where the auxiliary function zi h∈RT1is used in 2 iteration steps of the majorant minimization. We see that error control for the variable pdominates in this case and takes almost half of the computational time of entire step. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx 23 Table 11 Example 3. Comparison of the CPU time required to perform several stages of the numerical calculations w.r.t. the selected time steps, h=1 64 . ncomputation of p, s (%) computation of u, s (%) error control for p, s (%) error control for u, s (%) (a) N=100 (a-1) 2 iterations for Mh pminimization, zi h∈RT1 20 96.05 (8.46%) 305.87 (26.94%) 544.04 (47.91%) 189.56 (16.69%) 40 100.00 (8.38%) 305.29 (25.58%) 539.27 (45.18%) 249.01 (20.86%) 80 103.62 (8.60%) 308.15 (25.56%) 561.49 (46.58%) 232.13 (19.26%) (a-2) 0 it. for Mh pminimization, zi h∈RT0 20 35.58 (21.62%) 118.02 (71.69%) 2.50 (1.52%) 8.51 (5.17%) 40 37.70 (22.35%) 123.79 (73.38%) 2.10 (1.25%) 5.10 (3.02%) 80 31.03 (21.24%) 106.15 (72.66%) 2.51 (1.72%) 6.40 (4.38%) (b) N=10 (b-1) 2 iteration for Mh pminimization, zi h∈RT1 2 29.96 (12.45%) 99.34 (41.27%) 44.11 (18.33%) 67.30 (27.96%) 4 35.43 (14.19%) 113.12 (45.31%) 44.51 (17.83%) 56.58 (22.66%) 8 48.91 (13.21%) 163.18 (44.09%) 59.82 (16.16%) 98.22 (26.54%) (b-2) 0 iteration for Mh pminimization, zi h∈RT0 2 43.41 (21.09%) 148.86 (72.30%) 2.92 (1.42%) 10.69 (5.20%) 4 36.39 (21.29%) 121.94 (71.35%) 2.36 (1.38%) 10.21 (5.98%) 8 40.92 (21.28%) 138.38 (71.96%) 2.65 (1.38%) 10.33 (5.37%) After that, we illustrate the CPU costs for the function zi h∈RT0without optimization procedure. Columns four and five emphasize how inexpensive a posteriori error control can be when one only needs to calculate Mh pand Mh u. 8. Conclusions and future work We analyze semi-discrete approximations of the Biot poroelastic problem and deduce guaranteed and fully computable bounds of corresponding errors. The derivation combines the estimates for contraction mappings and the functional a posteriori error majorants for elliptic problems. The obtained error bound is fully computable and independent on the discretization techniques used for the variational formulation of the Biot problem as soon as the reproduced approximations belong to admissible functional spaces. Moreover, obtained error functional does not depend on any mesh discretization constants and only contains global Poincare-type constants characterizing considered geometry. Numerical results presented above provided the evidence of the quantitative efficiency of the majorant when it comes to the error indication. Generally, automation of the mesh adaptation procedure crucially depends on the quality of the error indicator used. It is a very important and complex topic that includes not only the question of how to reconstruct the mesh or which approximation space to use but also how to combine it with existing well-verified technologies (e.g., greedy marking, hp-refinement, etc.). Thorough investigation of these questions in the context of the error estimates introduced in this paper is an important task for future research. CRediT authorship contribution statement Kundan Kumar: Conceptualization, Formal analysis, Funding acquisition, Methodology, Writing - review & editing. Svetlana Kyas: Conceptualization, Formal analysis, Methodology, Software, Investigation, Validation, Writing - original draft, Writing - review & editing. Jan Martin Nordbotten: Conceptualization, Formal analysis, Methodology, Writing - review & editing. Sergey Repin: Conceptualization, Formal analysis, Methodology, Writing - review & editing. Acknowledgments The work is funded by the SIU Grant CPRU-2015/10040. The second author wants to acknowledge the Werner Siemens Foundation (Werner Siemens-Stiftung) for its support of the Geothermal Energy and Geofluids group at ETH Zurich, Switzerland. KK and JMN would like to acknowledge Norwegian Research Council project Toppforsk 250223 for funding. Pleasecitethisarticleas:K. Kumar,S. Kyas,J.M. Nordbottenetal.,Guaranteedandcomputableerrorboundsforapproximationsconstructedbyaniterative decoupling of the Biot problem, Computers and Mathematics with Applications (2020), https://doi.org/10.1016/j.camwa.2020.05.005. 24 K. Kumar, S. Kyas, J.M. Nordbotten et al. / Computers and Mathematics with Applications xxx (xxxx) xxx Appendix Notation and definitions of spaces. We use the standard Lebesgue space of square-measurable functions L2(Ω) equipped with the norm ∥v∥Ω:= ∥v∥L2(Ω):= (v, v)1/2 Ωfor all u, v ∈L2(Ω). Let Md×ddenote the space of real d-dimensional tensors. The products of vector-valued v,w∈Rdand tensor-valued functions τ,σ∈Md×dare defined by the relations (v,w)Ω:= ∫Ω v·wdxand (τ,σ)Ω:= ∫Ω τ:σdx, where v·w:= viwiand τ:σ:= τij σij, respectively. Next, A(x)∈Md×d,x∈Ωdenotes a symmetric uniformly positive defined matrix that satisfies 0 < λ ≤λ(x)≤λ≤ +∞,λ, λ ∈Rwith uniformly bounded eigenvalues λ(x). Then, for the product (u,v)A:= (Au,v), we have λ∥v∥≤∥v∥A≤λ∥v∥and (u,v)≤ ∥·∥A∥·∥A−1,∀u,v∈ [L2(Ω)]d. We use the standard notation for the Sobolev space of vector-valued functions having square-summable derivatives H1(Ω):= {v∈L2(Ω)| ∇v∈ [L2(Ω)]d}, equipped with the norm ∥v∥H1(Ω):= (∥v∥2 Ω+ |v|2 Ω)1/2. Also, we use a semi-norm |v|Ω:= |v|H1(Ω):= ∥∇ v∥Ω, and for the vector-valued functions with square-summable divergence introduce the Hilbert space H(Ω,div) := {v∈ [L2(Ω)]d|divv∈L2(Ω)}, endowed with the norm ∥v∥2 H(Ω,div) := ∥v∥2 [L2(Ω)]d+∥divv∥2 L2(Ω). Let Σbe a part of the boundary such that measd−1Σ>0 (in particular, it may coincide with ∂Ω). For the functions in H1 0,Σ(Ω):= {v∈H1(Ω)|v|Σ=0}, the Friedrichs-type inequality reads: ∥v∥Ω≤CF Γ|v|Ω,∀v∈H1 0,Σ(Ω).(7.2) The corresponding trace operator γ:H1(Ω)→H 1 2(Ω) is bounded and satisfies the estimate v|Σ:= γ v, ∥v∥Σ≤Ctr ΣΩ ∥v∥H1(Ω),∀v∈H1(Ω),(7.3) where ∥v∥Σis the norm of L2(Σ). CKdenotes the constant in the Korn inequality ∥w∥[H1(Ω)]d≤CK∥ε(w)∥[L2(Ω)]d×d,∀w∈ [H1(Ω)]d,(7.4) Also, we use the inequality ∥div w∥=∥tr ε(w)∥=∥I:ε(w)∥ ≤ √d∥ε(w)∥,∀w∈ [H1(Ω)]d,(7.5) where I∈Md×dis the unit tensor of Md×d, and ε(w)∈Md×ddenotes the symmetric part of ∇w. Next, let Q:= Ω×(0,T) denote a space–time cylinder (with given time-interval (0,T), 0 <T<+∞), and let Σ=∂Ω×(0,T) be a lateral surface of Q, whereas Σ0:= ∂Ω×{0}and ΣT:= ∂Ω×{T}define the bottom and the top part of the mantel (so that ∂Q=Σ∪Σ0∪ΣT). Consider functions defined in (0,T) with values in a functional space X (cf. [71–73]). Let ∥·∥Xdenote the norm in X, then for r=2, we define the Bochner space L2(0,T;X):= {fmeasurable in [0,T]⏐⏐⏐∫T 0∥f(t)∥2 Xdt<∞}, and respective norm ∥f∥L2(0,T;X):= (∫T 0∥f(t)∥2 Xdt)1/2. It is a Hilbert space if Xis a Hilbert space. Throughout the paper, we also use the spaces H1(0,T;X):= {f∈L2(0,T;X)|∂tf∈L2(0,T;X)}(7.6) equipped with norm ∥u∥H1(0,T;X):= (∫T 0(∥∂tf(t)∥2 X+∥f(t)∥2 X)dt)1/2. We assume that This a regular mesh satisfying angle condition defined on Ω. Then, the corresponding discretization spaces with the Lagrangian finite elements of order 0 or 1 are defined as P0:= {vh∈L2(Ω)|∀T∈Th, vh|T∈P0},P1:= {vh∈H1(Ω)|∀T∈Th, vh|T∈P1},(7.7) where Pkdenotes the space of polynomials of the order k∈N∪0. The Raviart–Thomas elements of the lowest and first order are denoted by RT0:= {yh∈H(div,Ω): ∀T∈Th,yh|T=a+bx,a∈Rd,b∈R}, RT1:= {yh∈H(div,Ω): ∀T∈Th,yh(x)|T=q(x)+xr(x),q∈ [P1]d,r∈P1}, respectively. Finally, the table below presents notation used in the paper for the physical quantities (see Table A.12).