A posteriori error estimates for hierarchical mixed-dimensional elliptic equations
Jhabriel Varela, Elyes Ahmed, Eirik Keilegavlen, Jan Martin Nordbotten, Florin Adrian Radu
AA posteriori error estimates for hierarchicalmixed-dimensional elliptic equations
Jhabriel Varela ∗ , † Elyes Ahmed ‡ Eirik Keilegavlen † Jan Martin Nordbotten † Florin Adrian Radu † January 22, 2021
Abstract
In this paper, we derive a posteriori error estimates for mixed-dimensionalelliptic equations exhibiting a hierarchical structure. Exploiting the exteriorcalculus perspective of such equations, we introduce mixed-dimensional vari-ables and operators, which, together with careful construction of the functionalspaces, allow us to recast the set of partial differential equations as a regularlinear elliptic problem structure-wise. We therefrom apply the well-establishedtheory of functional a posteriori error estimates to our model to derive guar-anteed abstract as well as fully computable upper bounds. Our estimators aretested using three different families of locally-mass conservative methods on syn-thetic problems and verification benchmarks of flow in fractured porous media.The numerical results support our theoretical findings while showcasing satis-factory effectivity indices.
Keywords : Mixed-dimensional geometry, A posteriori error estimates, Frac-tured porous media.
Mixed-dimensional partial differential equations (mD-PDEs) arise when partial differ-ential equations interact on domains of different topological dimensions [36]. Proto-typical examples include models of thin inclusions in elastic materials [6, 17, 12], bloodflow in human vasculature [20, 31, 27], root water uptake systems [30], and flow infractured porous media [5, 24, 2]. The latter example has an appealing mathematicalstructure, in that the model equations allow for a hierarchical representation whereeach domain (matrix, fractures, and their intersections) only has direct interactionswith domains of topological dimension one higher or one lower [14]. Such hierarchicalmD-PDEs are the topic of the current paper.mD-PDEs are intrinsically linked to the underlying geometric representation,which, in a certain sense, generalizes the usual notion of the domain. One can thendefine sets of suitable functions (and function spaces) on this geometry, and thesesets are then naturally interpreted as mixed-dimensional functions. Exploiting this ∗ This work was supported by the Norwegian Academy of Science and Letters through the VISTAproject 6371. † Department of Mathematics, University of Bergen, P.O. Box 7800, N-5020 Bergen, Norway.Corresponding author e-mail: [email protected]. ‡ SINTEF Digital, Mathematics and Cybernetics, P.O. Box 124 Blindern, N-0314 Oslo, Norway. a r X i v : . [ m a t h . NA ] J a n oncept, one can generalize the standard differential operators to mappings betweenmixed-dimensional functions and obtain a mixed calculus.The fact that this mixed-dimensional calculus inherits standard properties of cal-culus, particularly integration-by-parts (relative to suitable inner products), a deRham complex structure, and a Poincar´e-Friedrichs inequality, was recently estab-lished using the language of exterior calculus on differential forms [13]. These toolsare the key to unlocking abstract a posteriori error bounds for mD-PDEs, followinga functional approach [45].This paper exploits the exterior calculus perspective of mD-PDEs to obtain aposteriori error guaranteed upper bounds for the mixed-dimensional Laplace equation.This corresponds to (in some cases, a generalization of) standard mixed-dimensionalformulations used for flow in fractured porous media. With this example in mind, wemake the abstract error bounds concrete and directly interpret the error bounds interms of classical function spaces. As is customary for functional-based a posteriori error bounds, the derivation is agnostic concerning the choice of how the approximatesolution is obtained. The error bounds are validated against analytical solutionsusing approximate solutions obtained by various numerical methods and display verysatisfactory effectivity indices.A wide variety of a posteriori error techniques are currently available for fixed-dimensional elliptic problems; e.g.: residual type a posteriori error estimates [7, 49],methods based on post-processing [54, 7], goal-oriented type estimates [40], and thefunctional type [43, 44, 46]. Nevertheless, existing a posteriori error bounds for mixed-dimensional models are available only for specific cases in the context of mortar meth-ods [53, 9, 52] and fractured porous media [16, 33, 26], with much less geometricgenerality than what we present herein. Thus, for practical problems, a posteriori error bounds have until now essentially not been available. The current work recti-fies this situation. We illustrate our error bounds’ applicability on recent communitybenchmark problems [23, 10] while considering three different families of locally-massconservative discretization methods.The paper is structured as follows: In section 2, we begin by introducing themixed-dimensional geometric decomposition, the model problem, an alternative com-pact notation, and the functional spaces in a continuous setting. We finish the sectionproviding the abstract a posteriori error bounds. Section 3 introduces the grid par-titions, functional spaces, and the model problem in a finite sense. In section 4,we derive the computable bounds and introduce the scaled version of the majoranttogether with the identification of local errors. The section is finished by providingconcrete applications of our estimators to locally-mass conservative methods. Sections5 and 6 deal, respectively, with validations and practical applications of the derivederror bounds. Finally, in section 7, we present our conclusions. In this section, we start with the description of the mixed-dimensional geometry andthe underlying physical model. We then introduce the mathematical quantities andfunction spaces used to characterize the weak solution. Therefrom, we infer abstract aposteriori error estimates featuring guaranteed upper bounds, and a global efficiencyresult, in the context of mixed-dimensional equations.2 Ω Ω Ω Γ + Γ − Γ Γ Figure 1: Mixed-dimensional geometric entities. Left: A one-dimensional fracture Ω with zero-dimensional tips Ω and Ω is embedded in a matrix Ω . Right: Interfacescoupling subdomains one dimension apart. We use the positive and negative signs todistinguish the different sides of Γ . We start by providing the mixed-dimensional geometric decomposition following [14,38, 13]. Consider an n − dimensional domain Ω decomposed into subdomains of dif-ferent dimensionalities (see Figure 1). The ambient dimension n can be either 2 or 3.The subdomains of dimension n − i and a dimension d i = d ( i ). That is, Ω d i i refers to the domain number i ∈ I with dimensionality d i = { , , . . . , n } . We consider a total of m subdomains, i.e., I = { , , . . . , m } , and referto the union of all subdomains as Ω = (cid:91) i ∈ I Ω d i i . (2.1)Intersections between subdomains of equal dimensionality are considered to be ( d − d –dimensional subdomains are disjointed.We allow subdomains one dimension apart to be connected via interfaces Γ d i ij , whichare defined by being adjacent to a lower-dimensional subdomain Ω d i i and a higher-dimensional subdomain Ω d j j , where d j = d i + 1. We will consistently use the firstsub-index to refer to the lower-dimensional subdomain and the second to the higher-dimensional subdomain. Occasionally, we will need to identify the sides of an interface.To avoid introducing extra notation, we will identify the interface’s sides as positiveand negative, i.e., Γ d i ij = Γ d i ij + ∪ Γ d i ij − (see Figure 1). Since there is no room forambiguity in what follows, we will drop the superscript denoting the dimensionalityof subdomains, interfaces, and mixed-dimensional variables in general.To establish the connection between subdomains and interfaces, we introduce twosets of global indices, namely ˆ S i and ˇ S i . The set ˆ S i contains the global indices ofthe higher-dimensional neighboring subdomains of co-dimensionality one w.r.t. Ω i .Conversely, the set ˇ S i contains the indices of the lower-dimensional neighboring sub-domains w.r.t. Ω i . For convenience, we define the union of all interfaces asΓ = (cid:91) i ∈ I (cid:91) j ∈ ˇ S i Γ ji , (2.2)3 ij ∂ i Ω j Ω i Ξ jij Ξ iij Π jij Π iij Figure 2: Generic coupling between a lower-dimensional subdomain Ω i and a higher-dimensional subdomain Ω j through a common interface Γ ij . The projector Ξ mapsfluxes from the interface to the subdomains, whereas the projector Π maps fluxesfrom subdomains to the interface. The part of the boundary of Ω j that geometricallycoincides with Γ ij is denoted by ∂ i Ω j .and identify the boundary of the domain Ω as ∂ Ω = (cid:91) i ∈ I ∂ Ω i \ Γ . (2.3)We will refer to the disjointed union of Ω and Γ as F , i.e., F = Ω (cid:116) Γ, and ∂ F = ∂ Ω.Now, considering a single interface Γ ij as shown in Figure 2, we define the normalvector n j on ∂ i Ω j to be directed from Ω j towards Ω i , where ∂ i Ω j represents the partof ∂ Ω j geometrically coinciding with Γ ij . Note that, in practice, the three geometricalobjects, namely ∂ i Ω j , Γ ij , and Ω i , coincide spatially. However, the distinction betweenthem is useful to consistently keep track of the transferring of variables [14].We can now focus on the decomposition of the boundaries of the domain. Weassume that the boundary of Ω is partitioned such that ∂ Ω = ∂ Ω D ∪ ∂ Ω N , with ∂ Ω D ∩ ∂ Ω N = ∅ , and ∂ Ω D with positive measure. Moreover, we assume that each n –dimensional subdomain Ω i of the surroundings has a non-empty Dirichlet boundarycondition, i.e., | ∂ Ω i ∩ ∂ Ω D | >
0, for all i with d i = n . We also introduce the portionof the boundaries of the subdomains corresponding to either a Dirichlet or a Neumannboundary condition. That is, for all i ∈ I , we have ∂ D Ω i = ∂ Ω i ∩ ∂ Ω D , ∂ N Ω i = ∂ Ω i ∩ ∂ Ω N , d i = { , . . . , n } . In what follows, whenever we explicitly restrict the topological scope of applica-bility of a mixed-dimensional variable or equation, i.e., using d i = {· , . . . , ·} , we willassume that such declaration holds for every i ∈ I . Before presenting the set of governing equations, we introduce the operators andprojectors in the mixed-dimensional context and the concept of mortar flux [14, 25].We will use the usual gradient and divergence operators ∇ ( · ) and ∇ · ( · ) acting onscalars and vectors, respectively. Note that if d i = { , . . . , n − } , the operators aredefined on their corresponding tangent spaces.Let us now introduce the projectors Ξ and Π as in [29]. The operators are usedto keep track of the mappings from the interfaces to the subdomains and vice versa(see Figure 2). The pair of subindices identifies the interface, whereas the superindex4dentifies the subdomain. We will use four types of projector operators:Ξ jij : Γ ij → ∂ i Ω j , Ξ iij : Γ ij → Ω i , (2.4)Π jij : ∂ i Ω j → Γ ij , Π iij : Ω i → Γ ij . (2.5)Naturally, in a continuous setting, we have that Ξ kij Π kij = I for k = { i, j } . However,this might not be the case in a discrete setting, where the operators become grid-dependent and possibly discretization-dependent, as discussed in section 3.3.To handle the inter-dimensional contributions to mass conservation, the amount ofmass entering through an interface is now introduced as an additional variable called mortar flux [14] : λ ij = Π jij tr u j · n j , on Γ ij , d i = { , . . . , n − } , (2.6)where u j is the velocity of the subdomain Ω j , n j is the normal vector on ∂ i Ω j (recallthat n j points from the higher- to the lower-dimensional subdomain), and tr is asuitable trace operator mapping from Ω j to ∂ i Ω j (see Figure 2).Since mass may enter the fracture from one side and continue tangentially throughthe fracture (creating a difference in normal fluxes) [39], we conveniently define anoperator capturing the jump across the interface Γ ij as in [14]:[[ · ]] i : Γ ij → Ω i , [[ λ ]] i = (cid:88) j ∈ ˆ S i Ξ iij λ ij , d i = { , . . . , n − } , (2.7)where λ ij represents the mortar flux associated with Γ ij . Note that for subdomainswith dimensionality equal to the ambient dimension, the jump in mortar fluxes iszero, i.e., [[ λ ]] i = 0 for all i with d i = n . We will now present the model problem that is valid for an arbitrary number offractures and intersections. This set of equations was previously proposed in [14, 38]and utilized in [29].Darcy’s law for an arbitrary subdomain reads u i = −K i ∇ p i , in Ω i , d i = { , . . . , n } , (2.8)where K i refers to the tangential component of the permeability tensor (satisfyingthe usual symmetry, positive definiteness, and uniformity requirements) and p i is thefluid pressure.The mass conservation equation is given by ∇ · u i − [[ λ ]] i = f i , in Ω i , d i = { , . . . , n } , (2.9)where the first term represents the divergence of the tangential velocity field, thesecond term takes into account the jump in mortar fluxes (2.7), and f i ∈ L (Ω i ) isan external source (or sink) term.We also need a constitutive relationship holding in the normal direction. Severalchoices are possible in this regard. We refer to [32, 37] for a thorough discussion inthis matter. To keep the model as general as possible, we assume that the pressurejump can be related to the mortar fluxes: λ ij = − κ ij (cid:16) Π iij p i − Π jij tr p j (cid:17) , on Γ ij , d i = { , . . . , n − } . (2.10)5ere, κ ij represents the normal effective permeability, and the terms in parenthesesdenote the jump in pressures perpendicular to Γ ij .The model is closed with proper internal and external boundary conditionstr u j · n j = Ξ iij λ ij , on ∂ i Ω j , d i = { , . . . , n − } , (2.11) p i = g Di , on ∂ D Ω i , d i = n, (2.12) u i · n i = g Ni , on ∂ N Ω i , d i = n, (2.13) u i · n i = 0 , on ∂ N Ω i , d i = { , . . . , n − } . (2.14)As commonly considered in practical applications, (2.14) imposes the so-calledimmersed tip boundary condition on fully embedded lower-dimensional subdomains[29, 25]. Equations (2.8) – (2.14) hold for all i ∈ I and represent the strong form ofthe model problem. This subsection will introduce a new notation that will allow us to recast the model,and perhaps most importantly, perform the analysis more compactly.First, let us define the mixed-dimensional pressure as the ordered collection ofpressures associated with each subdomain, namely p := [ p i ]. Similarly, we define themixed-dimensional flux as the doublet of tangential and normal fluxes u := [ u i , λ ij ].We will consistently use gothic fonts to refer to mixed-dimensional quantities. When-ever possible, we will use ‘mD’ instead of mixed-dimensional.Mimicking the procedure from [28, 25], we define the mD divergence D · ( · ) as D · u = D · [ u i , λ ij ] := [ f i ] . (2.15)Similarly, we define the mD gradient D ( · ) as D p = D [ p i ] := [ r i , r ij ] , (2.16)where r i = ∇ p i and r ij = Π iij p i − Π jij tr p j . (2.17)The material parameters are collected into an mD function K , such that K − [ u i , λ ij ] := (cid:2) K − i u i , κ − ij λ ij (cid:3) , (2.18)where we assume that K = K T . Moreover, we assume that K has the smallesteigenvalue λ K , such that ( K u , u ) L ( F ) ≥ λ K ( u , u ) L ( F ) , (2.19)where the inner product ( · , · ) L ( F ) will be formally defined in the next section.With the above definitions, we can recast our model problem as D · u = f , in Ω , (2.20a) u = − K D p , in Ω × Γ , (2.20b)where f = [ f i ] is the mD source term. Note that proper external boundary conditions(2.12)–(2.14) should complement the model.6 .5 Weak forms and functionals In this subsection, we will establish the weak forms of our model problem. To this aim,let us first introduce the functional spaces of interest. In a slight abuse of notation,we define the mD square-integrable space for the pressure as L ( F ) := L (Ω) = (cid:89) i ∈ I L (Ω i ) . (2.21)Similarly, we define the mD square-integrable space for the fluxes as L ( F ) := L (Ω) × L (Γ) = (cid:89) i ∈ I L (Ω i ) × (cid:89) i ∈ I (cid:89) j ∈ ˇ S i L (Γ ji ) , (2.22)where L (Ω i ) := L (cid:0) (Ω i ) d i (cid:1) . We will use the space of mD weakly-differentiablefunctions H ( F ) := (cid:8) q ∈ L ( F ) : D q ∈ L ( F ) (cid:9) , (2.23)and the space of mD weakly-differentiable functions with compact support on external Dirichlet boundaries˚ H ( F ) := { q ∈ H ( F ) : [ q i ] = 0 , ∀ i ∈ I | d i = n } . (2.24)To define weakly differentiable mD flux spaces, we first need to enhance the con-tinuity of the square integrable flux spaces. Therefore, we introduce the standardextension operator for Neumann boundary data from boundaries onto domains, as in[14]. We denote such operators as R ji : L (Γ ji ) → H (div , Ω i ), satisfying that for any λ ∈ L (Γ ji ) (tr R ji λ )( y ) = (cid:40) λ ( y ) y ∈ Γ ji . (2.25)The precise choice of R is not important; however, the natural choice based on thesolution of an auxiliary elliptic equation is reasonable [14]. This allows us to definean mD flux space as follows H ( D · , F ) := (cid:89) i ∈ I ˚ H (div , Ω i ) × (cid:89) i ∈ I (cid:89) j ∈ ˇ S i R ji L (Γ ji ) . (2.26)The space H ( D · , F ) has slightly enhanced regularity relative to a space composed of H (div , Ω i ) on each subdomain, as this latter case would not have traces in L (Γ ji ).An alternative construction of H ( D · , F ) can be obtained starting with spaces withstronger continuity and weakening their regularity based on the mD differential op-erators and integration-by-parts. Both constructions result in equivalent functionalspaces [13].Following [37], we introduce two types of mD inner products, namely one for themD pressures ( · , · ) L ( F ) and another for the mD fluxes ( · , · ) L ( F ) . The first one isgiven by the sum of the usual L -inner products of the subdomains( p , q ) L ( F ) = (cid:88) i ∈ I ( p i , q i ) L (Ω i ) , p , q ∈ L ( F ) . (2.27)The second inner product is given by the sum of the subdomains’ L -inner prod-ucts plus the traces of the lower-dimensional adjacent interfaces in the direction per-pendicular. That is, for any u , v ∈ L ( F )( u , v ) L ( F ) = (cid:88) i ∈ I ( u i , v i ) L (Ω i ) + (cid:88) j ∈ ˇ S i (cid:0) Π iji tr u i · n i , Π iji tr v i · n i (cid:1) L (Γ ji ) . (2.28)7t is easy to check that (2.27) and (2.28) naturally induce the following norms (cid:107) q (cid:107) L ( F ) = ( q , q ) / L ( F ) ∀ q ∈ L ( F ) , (2.29) (cid:107) v (cid:107) L ( F ) = ( v , v ) / L ( F ) ∀ v ∈ L ( F ) . (2.30)A key tool for our analysis is integration-by-parts in the context of mixed-dimensionaloperators. This has been formally proven in [13] and henceforth be stated here with-out proof. Lemma 2.1 (Integration-by-parts for mD operators) . Let v ∈ H ( D · , F ) and q ∈ ˚ H ( F ) arbitrary, then ( D · v , q ) L ( F ) = − ( v , D q ) L ( F ) . (2.31)The variational formulations are now straightforward following [14, 37, 38] butmore closely [13] using exterior calculus. The primal weak form consists of finding p ∈ ˚ H ( F ) such that ( K D p , D q ) L ( F ) = ( f , q ) L ( F ) ∀ q ∈ ˚ H ( F ) , (2.32)and the dual weak form consists of finding u ∈ H ( D · , F ) and p ∈ L ( F ) such that (cid:0) K − u , v (cid:1) L ( F ) − ( p , D · v ) L ( F ) = 0 ∀ v ∈ H ( D · , F ) , (2.33a)( D · u , q ) L ( F ) = ( f , q ) L ( F ) ∀ q ∈ L ( F ) . (2.33b)The next two functionals are necessary for the a posteriori analysis following [51].Let the symmetric bilinear form B ( · , · ) acting on mD pressures be defined as B ( p , q ) := ( K D p , D q ) L ( F ) , p , q ∈ ˚ H ( F ) , (2.34)whereas A ( · , · ) acting on mD fluxes be given by A ( u , v ) := (cid:0) u , K − v (cid:1) L ( F ) , u , v ∈ L ( F ) . (2.35)It is easy to see that the above mD bilinear forms induce the following permeability-weighted semi-norms on the space ˚ H ( F ): ||| q ||| := B ( q , q ) = (cid:13)(cid:13)(cid:13) K D q (cid:13)(cid:13)(cid:13) L ( F ) , ∀ q ∈ ˚ H ( F ) , (2.36)and the mD energy semi-norm ||| v ||| ∗ := A ( v , v ) = (cid:13)(cid:13)(cid:13) K − v (cid:13)(cid:13)(cid:13) L ( F ) , ∀ v ∈ L ( F ) . (2.37)The semi-norms (2.36) and (2.37) are related via ||| q ||| = ||| K D q ||| ∗ , ∀ q ∈ ˚ H ( F ) . (2.38) a posteriori error estimates Having available all the ingredients, in this subsection, we will obtain abstract boundsbased on functional a posteriori estimates [35, 45]. The following bounds are abstractbecause they are agnostic to any numerical method and no grid-dependent quantityneeds to be introduced. We summarize the results in the following theorem.8 heorem 2.2 (Abstract a posteriori error estimates) . Let p ∈ ˚ H ( F ) and u ∈ H ( D · , F ) be respectively the exact solutions to (2.32) and (2.33) , and let q ∈ ˚ H ( F ) and v ∈ H ( D · , F ) be arbitrary. Then, the error on the mD pressure measured in thenorm (2.36) is bounded such that ||| p − q ||| ≤ ||| v + K D q ||| ∗ + C F λ K (cid:107) f − D · v (cid:107) L ( F ) := M ( q , v ; f ) , (2.39) where C F is the Poincar´e-Friedrichs constant, and M ( q , v ; f ) is referred to as themajorant. Moreover, the error on the mD flux measured in the norm (2.37) is alsobounded by the majorant. That is, ||| u − v ||| ∗ ≤ M ( q , v ; f ) . (2.40) Proof.
The proof for the bounds on the pressure follows the classical case closely (see,e.g. [45]) as generalized to the mD setting. Choose any q ∈ ˚ H ( F ) and measure thedifference between q and the exact solution p using the energy norm (2.36): ||| p − q ||| = B ( p − q , p − q ) = ( K D ( p − q ) , D ( p − q )) L ( F ) = ( − K D q , D ( p − q )) L ( F ) − ( u , D ( p − q )) L ( F ) = ( − K D q , D ( p − q )) L ( F ) + ( f , p − q ) L ( F ) . (2.41)Here, we used the definition of the bilinear form B ( · , · ), the constitutive law (2.20b),integration-by-parts (2.31), and (2.20a). Next, we use (2.31) to construct the identity( D · v , p − q ) L ( F ) + ( v , D ( p − q )) L ( F ) = 0 for an arbitrary v ∈ H ( D · , F ), and add itto (2.41), to get ||| p − q ||| = ( − ( v + K D q ) , D ( p − q )) L ( F ) + ( f − D · v , p − q ) L ( F ) . (2.42)Now, we apply the Cauchy-Schwarz inequality to the first and second terms of(2.42) to get the following bound ||| p − q ||| ≤ ||| v + K D q ||| ∗ + C F λ F (cid:107) f − D · v (cid:107) L ( F ) , where we have applied the norm definition (2.37), and the Poincar´e-Friedrichs in-equality given by (cid:107) p − q (cid:107) L ( F ) ||| p − q ||| ≤ C F λ K , (2.43)which has been formally proved in [13]. This completes the proof for the abstractupper bound on the mD pressure.The proof for the bounds on the mD flux is given next. We remark that analternate proof based on a generalized abstract estimate (Theorem 6.1, [51]) can beused to obtain equivalent upper bounds after its generalization to the mD setting.We start by considering the flux − K D q obtained from the arbitrary potential q ,and its normed difference from the exact flux: τ = u + K D q ||| u + K D q ||| ∗ . We can now define a projected flux t ∗ as the solution to the one-dimensional mDellitpic equation: t ∗ = arg min t | τ ( t )= ± τ ||| v − t ||| ∗ , τ ( t ) := ( u − t ) / ||| u − t ||| ∗ . Using Galerkin orthogonality arguments, the projectedflux t ∗ satisfies the standard property A ( u − t ∗ , v − t ∗ ) = 0 . Clearly, using the triangle inequality we have that ||| u − v ||| ∗ ≤ ||| u − t ∗ ||| ∗ + ||| v − t ∗ ||| ∗ , (2.44)while again thanks to Galerkin orthogonality ||| u − t ∗ ||| ∗ = A ( u − t ∗ , u − t ∗ ) ||| u − t ∗ ||| ∗ = A ( v − u , ± τ ) . (2.45)Finally, using t ∗ = − K D q and (2.45), (2.44) can be bound it as follows ||| u − v ||| ∗ ≤ ||| v + K D q ||| ∗ + A ( v − u , ± τ )= ||| v + K D q ||| ∗ + (cid:12)(cid:12)(cid:12)(cid:12) A (cid:18) v − u , u + K D q ||| u + K D q ||| ∗ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) = ||| v + K D q ||| ∗ + (cid:12)(cid:12)(cid:12)(cid:12) A (cid:18) v − u , K D ( q − p ) ||| K D ( q − p ) ||| ∗ (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) = ||| v + K D q ||| ∗ + (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:18) D · ( v − u ) , q − p ||| K D ( q − p ) ||| ∗ (cid:19) L ( F ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ ||| v + K D q ||| ∗ + C F λ K (cid:107) f − D · v (cid:107) L ( F ) , where we used (2.20b), mD integration-by-parts (2.31), and lastly the Poincar´e-Friedrichs inequality (2.43).The upper bounds from Theorem 2.2 can be used to obtain a result for the effi-ciency of the combined error. Corollary 2.3 (Efficiency of the combined error) . The majorant M ( q , v ; f ) from The-orem 2.2 has a guaranteed efficiency index no worse than 3 for the combined error E p , u . That is, M ( q , v ; f ) ≤ E p , u ≤ M ( q , v ; f ) , (2.46) where E p , u := ||| p − q ||| + ||| u − v ||| ∗ + C F λ K (cid:107) f − D · v (cid:107) L ( F ) . (2.47) Proof.
Use the triangle inequality in the definition of the majorant (2.39), to get M ( q , v ; f ) ≤ ||| v − u ||| ∗ + ||| u + K D q ||| ∗ + C F λ K (cid:107) f − D · v (cid:107) L ( F ) ≤ ||| u − v ||| ∗ + ||| v + K D q ||| ∗ + 2 C F λ K (cid:107) f − D · v (cid:107) L ( F ) ≤ M ( q , v ; f ) + C F λ K (cid:107) f − D · v (cid:107) L ( F ) ≤ M ( q , v ; f ) . To obtain the final bound, we used (2.20b), the equivalence between norms (2.38),the bound (2.40), and once again, the definition of the majorant (2.39).Even though the above results are attractive from a theoretical standpoint, inpractice, we would like to have an upper bound that is fully computable after obtainingan approximated solution using a numerical method. This requires us to introducethe grid and the approximated spaces, which are the focus of the next section.10 Γ ij T Ω i T Ω j T ∂ i Ω j Figure 3: Schematic representation of the coupling between the higher-dimensionalgrid T Ω j , the mortar grid T Γ ij , and the lower-dimensional grid T Ω i . Left: Strictlymatching coupling between the grids. Right: Fully non-matching coupling betweenthe grids. This section discusses the discrete setting, introducing the grid partition, discretespaces, and the approximated solution. Finally, we will need in the a posteriori errorestimates an ˚ H ( F )-conforming pressure for which a pressure reconstruction procedureis proposed. We will denote by T Ω i the partition of a given subdomain Ω i and let T Ω = (cid:83) i ∈ I T Ω i represent the partition of all subdomains. Similarly, let T Γ ij be the partition of aninterface Γ ij (which will be referred to as the mortar grid), with T Γ = (cid:83) i ∈ I (cid:83) j ∈ ˇ S i T Γ ji representing the union of all mortar grids. Analogously to the continuous setting, wewill denote by T ∂ i Ω j ⊂ T Ω j the part of T Ω j spatially coinciding with T Γ ij . Then, themixed-dimensional grid is defined as T F = T Ω (cid:116) T Γ .The grids will be constructed such that the elements K of T Ω i are strictly non-overlapping simplices of dimension d K = { , . . . , n } (note that for K ∈ T Γ ij , d K = { , . . . , n − } ). In addition, we will need T V Ω i , denoting the set of elements K ∈ T Ω i sharing a common node V .Notice that we allow T Γ ij to be non-matching with respect to their neighboringgrids (see Figure 3). We will see that the estimates are still valid for this more genericcase, provided that the discrete counterparts of the projector operators from section3.2 are well-defined.Finally, we let h K denote the diameter of K , where h Ω i = max h K T Ω i , h Γ ij =max h K T Γ ij , h ∂ i Ω j = max h K T ∂ i Ω j , and h = max T F . In this subsection, we will introduce the finite element spaces necessary to write theapproximated problem. We first define a local space for the approximated pressures,boundary fluxes, and internal fluxes as Q i,h = (cid:8) q i,h ∈ L (Ω i ) : q i,h | K ∈ P k ( K ) , ∀ K ∈ T Ω i (cid:9) , d i = { , · · · , n } , (3.1)Λ ji,h = (cid:8) µ ji,h ∈ L (Γ ji ) : µ ji,h | K ∈ P k ( K ) , ∀ K ∈ T Γ ji (cid:9) , d i = { , . . . , n } , (3.2)and V i,h = (cid:8) v i,h ∈ L (Ω i ) : v i,h | K ∈ RTN n ( K ) , ∀ K ∈ T Ω i (cid:9) , d i = { , . . . , n } . (3.3)11ere, P k represents the space of polynomials of degree k ≤ n , for some n ∈ Z + , and RTN n represents the Raviart-Thomas(-Nedelec) space of vector functions of degree n [42, 34].With the grids and the discrete spaces formally defined, the discrete counterpartof the projector operators (2.4) are denote by˜Ξ jij : L (cid:0) T Γ ij (cid:1) → L (cid:0) T ∂ i Ω j (cid:1) , ˜Ξ iij : L (cid:0) T Γ ij (cid:1) → L ( T Ω i ) , ˜Π jij : L (cid:0) T ∂ i Ω j (cid:1) → L (cid:0) T Γ ij (cid:1) , ˜Π iij : L ( T Ω i ) → L (cid:0) T Γ ij (cid:1) . The construction of the projectors will depend on the type of variables being mapped.If projections of the lowest order are used, this implies that intensive variables (e.g.pressures) should be averaged, whereas extensive variables (e.g. fluxes) should beintegrated [29].The composite space for the mD pressure Q h ⊂ L ( F ) is defined such that Q h = (cid:89) i ∈ I Q i,h , d i = { , . . . , n } . (3.4)Similarly, we construct the space X h ⊂ H ( D · , F ) for the mD flux X h = (cid:89) i ∈ I ˚ H (div , Ω i ) × (cid:89) i ∈ I (cid:89) j ∈ ˇ S i R ji ˜Ξ iji Λ ji,h ∩ V i,h , d i = { , . . . , n } . (3.5)Finally, we require that the discrete function spaces satisfy Q i,h = ∇ · V i,h , d i = { , . . . , n } , (3.6)which can easily be obtained by choosing any mixed-finite element method pairs (see[11]). Note, however, that this condition restricts the application of our estimates tolocally-mass conservative methods.In what follows, we will use the following relations linking continuous and brokennorms. For any q = [ q i ] ∈ L ( F ), it holds (cid:107) q (cid:107) L ( F ) = (cid:88) i ∈ I ( q i , q i ) L (Ω i ) = (cid:88) i ∈ I (cid:88) K ∈T Ω i ( q i , q i ) L ( K ) , (3.7)and for any v = [ v i , µ ij ] ∈ L ( F ) we have (cid:107) v (cid:107) L ( F ) = (cid:88) i ∈ I (cid:88) K ∈T Ω i ( v i , v i ) L ( K ) + (cid:88) j ∈ ˇ S i (cid:88) K ∈T Γ ji ( µ ji , µ ji ) L ( K ) . (3.8) We are now in position to introduce the approximated problem. In order to satisfycondition (3.6), we are restricted to use the approximation to the dual problem (2.33).Thus, we are interested in finding the mD pair ( p h , u h ) ∈ Q h × X h , such that (cid:0) K − u h , v h (cid:1) L ( F ) − ( p h , D · v h ) L ( F ) = 0 ∀ v h ∈ X h , (3.9a)( D · u h , q h ) L ( F ) = ( f , q h ) L ( F ) ∀ q h ∈ Q h . (3.9b)12n this paper, we are restricted to numerical methods exhibiting locally-mass con-servative properties. In other words, the so-called conforming methods such as finiteelement methods (FEM) are beyond the scope of this work (the extension should bestraightforward after a proper reconstruction of the fluxes, e.g. since u h / ∈ H ( D · , F )).We will therefore make this requirement explicit in the following assumption. Assumption 3.1 (Local mass conservation property) . Let u h = [ u i,h , λ ij,h ] be asolution of the approximated dual problem (3.9) . Then, we say that a numericalmethod is locally-mass conservative if for all i ∈ I : ( ∇ · u i,h − [[ λ h ]] i , L ( K ) = ( f i , L ( K ) , ∀ K ∈ T Ω i . (3.10)We can see that, by solving (3.9) and provided that Assumption 3.1 is satis-fied, the mD approximated flux is readily available for use in Theorem 2.2, since u h ∈ H ( D · , F ). However, note that this is not the case for the approximated mDpressure, which is only piecewise constant on each subdomain. Consequently, weare required to reconstruct p h into a function exhibiting higher regularity than theoriginal approximation. The reconstructed pressure is defined as follows. Definition 3.2 (Reconstructed pressure) . Let ( p h , u h ) ∈ L ( F ) × H ( D · , F ) be asolution to the approximated dual problem (3.9) . Then, we will call reconstructed mDpressure to any function ˜ p h obtained from the duplet ( p h , u h ) , such that ˜ p h ∈ ˚ H ( F ) . (3.11)In practice, there is some freedom for obtaining ˜ p h . Arguably, the simplest op-tion is to perform an averaging of the piecewise-constant pressure approximation onlocal patches [43]. A more elaborated alternative consists of solving local Neumannproblems first to obtain a quadratic post-processed pressure (which is only locallyconforming) and then apply the so-called Oswald interpolator to obtain a conformingreconstruction [22, 51, 4, 3].In this paper, we opt for balancing the complexity of the techniques mentionedabove and propose a reconstruction method that uses the local flux field information(hence arguably retaining most of the higher-order reconstruction accuracy) whilemaintaining a piecewise linear reconstruction.To be precise, we use the flux u i,h to obtain a local pressure gradient, which isthen used to project p i,h | K ∈ P k ( K ) onto the grid nodes. If C denotes the centroidof K ∈ T Ω i , the cell-centered pressure gradient is given by ∇ φ i,h ( C ) = −K − i | K u i,h ( C ) . Then, the reconstructed pressure for an internal node V can be written as˜ p i,h ( V ) = (cid:80) K ∈T V Ω i | K | ( p i,h ( C ) + ∇ φ i,h ( C ) · δ ( V )) (cid:80) K ∈T V Ω i | K | , (3.12)where δ ( V ) = x ( V ) − x ( C ) is the vector defined from the node V to the cell-center C of K . Inside every element K ∈ T Ω i , ˜ p i,h is defined as the affine function satisfying(3.12) at the nodes of K . For an external Dirichlet boundary node, we use˜ p i,h ( V ) = g Di ( V ) , d i = n, (3.13)where we assume sufficient regularity at the Dirichlet boundaries. We shall apply thisprocedure to all Ω i , and then define the mixed-dimensional reconstructed pressure as˜ p h := [˜ p i,h ]. 13 A posteriori error estimates
In this section, we derive guaranteed a posteriori error estimates in standard notationstarting from the abstract estimates. From there, a scaled version of the majorant,together with the distinction of local errors, will be introduced. Finally, we concretelypresent how classical locally-mass conservative numerical methods fit our framework.
Here, we show how abstract a posteriori estimates from section 2.6 straightforwardly(1) deliver fully and locally computable error estimates and (2) distinguish the variouserror components in a local sense. To this aim, we rephrase the local Poincar´e-Friedrichs inequality (see, e.g. [41, 8]) in the mD setting. For any ϕ ∈ H ( K ): (cid:107) ϕ − π ( ϕ ) (cid:107) L ( K ) ≤ C P ,K h K (cid:107)∇ ϕ (cid:107) L ( K ) , K ∈ T Ω i , i ∈ I, (4.1)where π ( ϕ ) denotes the mean value of ϕ over K , and C P ,K is the local Poincar´econstant, which is equal to 1 /π whenever K is convex.Note that, since the majorant bounds the pressure and the flux errors, the followingtheorem provides a practical and computable quantity for both errors. Theorem 4.1 (Fully computable version of the majorant) . Let ˜ p h ∈ ˚ H ( F ) be thereconstructed mD pressure from Definition 3.2, u h ∈ H ( D · , F ) the approximated mDflux satisfying Assumption 3.1, and f ∈ L ( F ) the mD source term. Then, M (˜ p h , u h ; f ) ≤ M h , (4.2) where M h = (cid:88) i ∈ I (cid:88) K ∈T Ω i η ,K + (cid:88) j ∈ ˇ S i (cid:88) K ∈T Γ ji η ,K + (cid:88) i ∈ I (cid:88) K ∈T Ω i η ,K , (4.3) is the computable version of the majorant of the error. Here, η DF ,K is the diffusive-flux estimator: η DF ,K = (cid:13)(cid:13)(cid:13) K − i u i,h + K i ∇ ˜ p i,h (cid:13)(cid:13)(cid:13) L ( K ) , K ∈ T Ω i , (cid:13)(cid:13)(cid:13) κ − ij λ ij,h + κ ij (cid:16) ˜Π iij ˜ p i,h − ˜Π jij tr ˜ p j,h (cid:17)(cid:13)(cid:13)(cid:13) L ( K ) , K ∈ T Γ ij , (4.4) and η R ,K the residual estimator: η R ,K = C P ,K h K c K ,K (cid:107) f i − ∇ · u i,h + [[ λ h ]] i (cid:107) L ( K ) , K ∈ T Ω i , (4.5) where C P ,K is the local Poincar´e constant and c K ,K := (cid:13)(cid:13)(cid:13) K − i (cid:13)(cid:13)(cid:13) L ( K ) .Proof. In Theorem 2.2, choose q = ˜ p h and v = u h , and use (2.42) to obtain ||| p − ˜ p h ||| ≤ ( − ( u h + K D ˜ p h ) , D ( p − q )) L ( F ) ||| p − ˜ p h ||| + ( f − D · u h , s ) L ( F ) = T + T . s = ( p − ˜ p h ) / ||| p − ˜ p h ||| ∈ ˚ H ( F ). The term T canbe bounded as follows: T = (cid:16) − K − ( u h + K D ˜ p h ) , K D ( p − ˜ p h ) (cid:17) L ( F ) ||| p − ˜ p h ||| ≤ (cid:13)(cid:13)(cid:13) K − u h + K D ˜ p h (cid:13)(cid:13)(cid:13) L ( F ) = (cid:34) (cid:88) i ∈ I (cid:32) (cid:88) K ∈T Ω i (cid:13)(cid:13)(cid:13) K − i u i,h + K i ∇ ˜ p i,h (cid:13)(cid:13)(cid:13) L ( K ) + (cid:88) j ∈ ˇ S i (cid:88) K ∈T Γ ji (cid:13)(cid:13)(cid:13) κ − ji λ ji,h + κ ji (cid:16) ˜Π jji ˜ p j,h − ˜Π iji tr ˜ p i,h (cid:17)(cid:13)(cid:13)(cid:13) L ( K ) (cid:33)(cid:35) , where we used the Cauchy-Schwarz inequality, the definition of the norm (2.30), andthe equivalence between continuous and broken norms (3.8).Finally, using Assumption 3.1 and definition (3.7), T can be easily estimatedfollowing [51, 1]: T = (cid:88) i ∈ I ( f i − ∇ · u i,h + [[ λ h ]] i , s i ) L (Ω i ) = (cid:88) i ∈ I (cid:88) K ∈T Ω i ( f i − ∇ · u i,h + [[ λ h ]] i , s i ) L ( K ) = (cid:88) i ∈ I (cid:88) K ∈T Ω i ( f i − ∇ · u i,h + [[ λ h ]] i , s i − π ( s i )) L ( K ) ≤ (cid:88) i ∈ I (cid:88) K ∈T Ω i (cid:107) f i − ∇ · u i,h + [[ λ h ]] i (cid:107) L ( K ) (cid:107) s i − π ( s i ) (cid:107) L ( K ) ≤ (cid:88) i ∈ I (cid:88) K ∈T Ω i (cid:107) f i − ∇ · u i,h + [[ λ h ]] i (cid:107) L ( K ) C P ,K h K (cid:107)∇ s i (cid:107) L ( K ) ≤ (cid:88) i ∈ I (cid:88) K ∈T Ω i C P ,K h K (cid:13)(cid:13)(cid:13) K − i (cid:13)(cid:13)(cid:13) L ( K ) (cid:107) f i − ∇ · u i,h + [[ λ h ]] i (cid:107) L ( K ) ||| s ||| . Here, we first used the mD inner product for scalars (2.27), together with theCauchy-Schwarz inequality and the Poincar´e-Friedrichs inequality. The assertion ofthe theorem follows from definition (2.36) and the fact that ||| s ||| = 1. Remark 4.2 (Nature of the estimate from Theorem 4.1) . The majorant M h fromTheorem 4.1 has two contributions: the error incurred in satisfying the constitutiverelationships in the tangential and normal directions ( η DF ,K ) and the error in satis-fying the mass conservation equation ( η R ,K ). Since the latter depends on the externalsources, this estimator is usually refer to as the data oscillation error [51]. Note that,for locally-mass conservative schemes, the residual estimator η R ,K usually takes verysmall values, and if f ∈ L ( F ) and piecewise constant then η R ,K = 0 . Remark 4.3 (Sharpness of M h ) . The inequality M (˜ p h , u h ; f ) ≤ M h follows from theuse of the local Poincar´e constant C P ,K instead of the global Poincar´e constant C F from Theorem 2.2. Therefore, M h gives a bound that is sharper than M (˜ p h , u h ; f ) . .2 Scaling of the majorant and identification of local errors While Theorem 4.1 gives a fully computable estimate for the global fracture network,in real applications, it is common to encounter large material contrasts between sub-domains [39]. For example, when the fracture network exhibits highly conductivefeatures ( K , κ → ∞ ) or nearly impermeable barriers ( K , κ → M h usually over-penalizes the error in the former case andunder-penalizes the error in the latter. This issue may lead to erroneous physical in-terpretations, such as that the largest errors concentrate in regions of small pressurejumps. To tackle this problem, we introduce a permeability-weighted version of themajorant as defined next. The resulting error bounds are algebraically equivalent tothose derived in the preceding section Definition 4.4 (Scaled version of the majorant) . Let M h be the majorant from The-orem 4.1. Then, we will denote by M ε to the permeability-weighted majorant, suchthat M ε := ξ − M h , (4.6) where ξ := max i ∈ I ( β i , γ i ) , (4.7) with β i := (cid:18) max K ∈T Ω i (cid:107)K i (cid:107) L ( K ) (cid:19) and γ i := (cid:34) max j ∈ ˇ S i (cid:32) max K ∈T Γ ji (cid:107) κ ji (cid:107) L ( K ) (cid:33)(cid:35) . (4.8)Similarly, it is often desirable to measure the error not only for the global fracturenetwork but also for its individual components. For this purpose, let us introduce thedistinction between the contribution to the global error associated with subdomainsand interfaces. Definition 4.5 (Scaled local errors) . Let ε Ω i and ε Γ ij be respectively the subdomainand interface errors given by ε Ω i := (cid:0) ε , Ω i + ε , Ω i (cid:1) , ε Γ ij := ε DF , Γ ij , (4.9) where ε , Ω i := (cid:88) K ∈T Ω i (cid:107)K i (cid:107) − L ( K ) η ,K , ε , Γ ij := (cid:88) K ∈T Γ ij (cid:107) κ ij (cid:107) − L ( K ) η ,K , (4.10) are respectively the subdomain and interface diffusive errors, and ε , Ω i := (cid:88) K ∈T Ω i η ,K , (4.11) is the subdomain residual error. It is easy to check that ε DF , Ω i , ε DF , Γ ij , and ε R , Ω i are related to M ε via M ε ≤ (cid:88) i ∈ I ε , Ω i + (cid:88) j ∈ ˇ S i ε , Γ ji + ξ − (cid:32)(cid:88) i ∈ I ε , Ω i (cid:33) . .3 Application to numerical methods In this subsection, we apply our estimates to three different families of classical locally-mass conservative methods that satisfy Assumption 3.1. These are; mixed finiteelement methods (MFEM), mixed virtual finite element methods (MVEM), and cell-centered finite volume methods (CCFVM).
MFEM aims at approximating the dual week problem (3.9) [11, 21]. For simplicity,we use only MFEM of the lowest order, e.g., Raviart-Thomas elements for fluxes andpiecewise constant approximations for the pressure and mortar fluxes. However, weremark that higher-order approximations fall into our framework without any extraconsiderations.For the lowest-order approximation, we use P ( K ) to build the spaces Q i,h andΛ ji,h from (3.1) and (3.2), and RTN ( K ) to construct the space V i,h from (3.3). Aftersolving (3.9) and reconstructing p h , the duplet (˜ p h , u h ) ∈ ˚ H ( F ) × H ( D · , F ) is readilyavailable for use in Theorem 4.1. Utilization of MVEM techniques has been steadily increasing due to their flexibility forapplication on general polygonal and polyhedral grids. For an excellent introductionto the method, we refer to [15]. Herein, we restrict our application to virtual elementsof the lowest order on simplicial elements.In this particular case, the tangential fluxes u i,h | K ∈ RTN ( K ) [19] and the normalfluxes λ ij,h | K ∈ P . Hence, we have that u h ∈ X h ⊂ H ( D · , F ). The method iscomplemented by a P ( K ) approximation of the pressures p h | K [25].After performing the pressure reconstruction, the duplet (˜ p h , u h ) ∈ H ( D · , F ) × ˚ H ( F ) can be applied to our framework. We must remark that, even though theapproximation spaces coincide for MFEM and MVEM in the case of simplicial mesheswith strictly non-overlapping elements, the actual discretization is different. For lowest order CCFVM in a mixed-dimensional context, pressures and mortar fluxesare approximated using P polynomials. However, unlike MFEM, the fluxes are notimmediately available for use in Theorem 4.1. In turn, we count with finite approxi-mations of the edge (face in 3d) fluxes satisfying the following local mass-conservativeproperty (cid:88) e ∈E K F K,e = ( f i , L ( K ) , ∀ K ∈ T Ω i , d i = { , . . . , n } , (4.12)where F K,e is the approximated flux trough the edge e of the element K , and E K isthe set of edges of the element K . Note that F K,e represents both standard fluxesand, for faces on the boundary to lower-dimensional objects, projections of mortarfluxes.Now the fluxes F K,e can be extended using
RTN on each element for each sub-domain [18, 50]. Finally, we can construct u h by collecting all the tangential andnormal fluxes, i.e., u h = [ u i,h , λ ij,h ]. After reconstructing p h , the duplet (˜ p h , u h ) ∈ ˚ H ( F ) × H ( D · , F ) falls naturally into the framework of Theorem 4.1.17 . . . . . Ω Ω . . . Ω . . ⌦
This section demonstrates our estimators’ capability by performing an efficiency anal-ysis using four different numerical methods, namely RT0-P0, MVEM, TPFA, andMPFA. The numerical examples were implemented in the open Python-based soft-ware
PorePy [29], using the extension package mdestimates [48], which includes thescripts of the numerical examples considered herein.We will use the methods mentioned above to approximate the solution of a 1D/2Dproblem with non-matching grids (section 5.1) and a 2D/3D problem with matchinggrids (section 5.2). The geometric configurations and fractures’ coordinates for bothproblems are shown in Figure 4.We will now proceed to describe the validation cases. For this purpose, let usdenote the fracture as Ω , the matrix as Ω , and the interfaces as Γ − and Γ + ,with Γ = Γ − ∪ Γ + .We start by assuming the existence of an exact pressure solution p in Ω , whichis given by the shortest distance from any point in Ω to Ω . This results in a piecewisequadratic function with three distinct regions for the 1D/2D problem (Equation (8.1))and nine regions for the 2D/3D problem (Equation (8.4)).Next, we set K = 1 and apply (2.8) to obtain the exact tangential velocity in Ω (see (8.2) and (8.5)). In particular, we note that u ( x ) = (cid:40) [1 , , . . . , T , < x < , ≤ x i ≤ ∀ i = { , . . . , d } , [ − , , . . . , T , < x < , ≤ x i ≤ ∀ i = { , . . . , d } , which implies that λ − = λ + = 1 after (2.6). Moreover, by setting κ − = κ + = 1,the pressure in the fracture is constrained by (2.10) such that p = −
1, meaning that u = [0 , . . . , T .The exact external source terms f ( x ) and f are now straightforward to obtainafter applying (2.9) (see (8.3) and (8.6)). The models are closed by imposing Dirichletboundary conditions satisfying the analytical pressure at the boundaries of Ω , andno-flux at the boundaries of Ω .We can now use the derived source terms and boundary conditions as input datafor approximating the solutions to the problems. To quantify the performance of The exact expressions of all relevant quantities involved in the validations can be found in theAppendix. ζ † ε Ω2 ε Ω1 ε Γ12 M h ||| p − ˜ p h ||| ||| u − u h ||| ∗ I eff p I eff u I eff p , u TPFA ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ ζ † The triplet ζ = ( h ∂ , h Γ12 , h Ω1 ) characterizes the non-matching coupling. The combination of meshsizes are: ζ = (0 . , . , . ζ = (0 . , . , . ζ = (0 . , . , . ζ = (0 . , . , . ζ = (0 . , . , . the numerical methods, we will measure the efficiency of the estimates using twoeffectivity indices, namely one for the pressure and one for the flux: I eff p = M h ||| p − ˜ p h ||| ≥ , and I eff u = M h ||| u − u h ||| ∗ ≥ . (5.1)In addition, we define a measure for the efficiency of the combined error as I eff p , u :=3 M (˜ p h , u h ; f ) / E p , u , which thanks to Corollary 2.3, it is guaranteed to lie between thefollowing bounds 1 ≤ I eff p , u ≤ . (5.2) For our first validation, we consider the 1D/2D configuration, as shown in the leftFigure 4. We are interested in the local and global estimates as well as in the effectivityindices I eff p , I eff u , and I eff p , u . Note that, in this particular case, M h = M (cid:15) since K = 1.For this problem, we assume that there is a non-matching coupling between T ∂ Ω , T Γ , and T Ω as depicted schematically in the right Figure 3. To characterize suchcoupling, we introduce the triplet ζ := ( h ∂ Ω , h Γ , h Ω ). We consider five levels ofsuccessively refined combinations of mesh sizes ζ = { ζ , . . . , ζ } as given in Table 1,where we also show the results.We first notice that all numerical methods show a clear local and global conver-gence tendency. The exception is TPFA, which exhibits issues for local convergence onthe matrix and fracture subdomains. This behavior goes in agreement with previouslyreported problems with TPFA [47]. We also note that the errors are quantitativelysimilar for MPFA, RT0-P0, and MVEM. This similarity has also been previouslyencountered [29].The columns I eff p and I eff u show the effectivity indices for the pressure and velocity.With exception to TPFA, the methods show an effectivity index of around 1 .
20 forthe pressure and 1 .
50 for the velocity. In the case of TPFA, the effectivity indices19able 2: Efficiency analysis for the 2D/3D validation problem. h † ε Ω2 ε Ω1 ε Γ12 M h ||| p − ˜ p h ||| ||| u − u h ||| ∗ I eff p I eff u I eff p , u TPFA h h h h h h h h h h h h h h h h h h h h † The mesh sizes are: h = 0 . h = 0 . h = 0 . h = 0 . h = 0 . show slightly different and more erratic values. In the last column, we see that I eff p , u related to the combined error remains nearly constant and lies within the expectedrange.Finally, we must remark to the reader that, as tempting as it might be to quantita-tively compare between the local errors, it is only fair to compare errors of subdomainsand interfaces of the same dimensionality. For our next numerical validation, we employ the 2D/3D configuration from the rightFigure 4. On this occasion, we use strictly matching grids for simplicity and againinvestigate five refinement levels, with the maximum mesh size given by h . The resultsare shown in Table 2.We can see that (up to the actual numerical values) the estimators give the sameresults as in section 5.1 in terms of local and global convergence. Again, TPFAstruggles to converge in the matrix and the fracture subdomains.For MPFA, RT0-P0, and MVEM, the effectivity indices for the pressure are closeto 1 .
1, and to 2 for the velocity. In the case of TPFA, the effectivity indices for thepressure are slightly larger than the other methods, whereas the velocity indices arecomparatively smaller. Once again, I eff p , u remains virtually constant and approximatelyequal to 2 . In this section, we apply our estimators to previous benchmark problems that havesolved the incompressible flow equations in fractured porous media.20 =
Mesh † ε Ω2 ε Ω1 , C ε Ω1 , B ε Γ1 , C ε Γ1 , B ε Γ0 M ε TPFA C 7.52e-1 3.75e-1 3.12e+0 2.28e-1 1.13e-1 5.85e-1 7.20e-2I 6.09e-1 2.74e-1 2.81e+0 1.38e-1 7.07e-2 4.51e-1 4.35e-2F 4.45e-1 1.75e-1 2.94e+0 7.11e-2 3.60e-2 3.22e-1 2.29e-2MPFA C 7.39e-1 3.83e-1 2.93e+0 2.16e-1 1.10e-1 5.55e-1 6.89e-2I 5.98e-1 2.78e-1 2.80e+0 1.34e-1 6.89e-2 4.43e-1 4.27e-2F 4.33e-1 1.78e-1 2.73e+0 6.71e-2 3.49e-2 3.00e-1 2.17e-2RT0-P0 C 7.39e-1 3.76e-1 4.33e+0 2.23e-1 1.12e-1 5.32e-1 7.11e-2I 5.95e-1 2.74e-1 4.14e+0 1.38e-1 7.12e-2 4.42e-1 4.41e-2F 4.30e-1 1.75e-1 4.10e+0 6.93e-2 3.60e-2 3.01e-1 2.24e-2MVEM C 7.29e-1 4.07e-1 2.37e+0 2.22e-1 1.14e-1 6.00e-1 7.03e-2I 5.91e-1 2.90e-1 2.35e+0 1.38e-1 7.22e-2 4.72e-1 4.39e-2F 4.28e-1 1.86e-1 2.31e+0 6.93e-2 3.62e-2 3.08e-1 2.23e-2 † C = Coarse, I = Intermediate, and F = Fine.
In this numerical experiment, we apply our estimators to the benchmark case 3b from[23]. As shown in Figure 5, the domain consists of ten (partially intersecting) fracturesembedded in a unit square matrix. The exact fracture coordinates can be found inAppendix C of [23]. The green lines represent conductive fractures ( K = 10 and κ = 10 ), whereas the red lines represent blocking fractures ( K = 10 − and κ = 1).The matrix permeability is set to one. A linear pressure drop is imposed from theleft ( p = 4) to the right ( p = 1), whereas no flux is prescribed at the top and bottomof the domains (see the right Figure 5).The benchmark establishes three refinement levels; coarse, intermediate, and fine,with approximately 1500, 4200, and 16000 two-dimensional cells. We show the localerrors for the matrix for two refinement levels obtained with MPFA in Figure 5.In Table 3, we show the local and global errors for the three refinement levels.To avoid numbering domains and interfaces, we refer to the matrix error as ε Ω , andgroup the fracture and interface errors by conductive and blocking. For example, ε Ω , C refers to the sum of the errors of (one-dimensional) conductive fractures.We note that all numerical methods show global convergence, and all except TPFAfor blocking fractures show local convergence. Note, in particular, that the conver-gence of blocking fractures is considerably slower than the conductive ones. This isa well-known limitation of the model that struggles to accurately represent nearly21 D F , K ( D )
Mesh † ε Ω3 ε Ω2 ε Ω1 ε Γ2 ε Γ1 ε Γ0 M ε TPFA C 9.43e-1 1.13e-5 2.34e-5 9.54e-2 8.09e-6 8.76e-6 4.53e-2I 4.82e-1 8.74e-6 1.20e-5 1.23e-2 2.71e-6 3.91e-6 5.67e-3F 4.12e-1 7.88e-6 9.14e-6 3.56e-3 1.42e-6 2.41e-6 1.57e-3MPFA C 6.46e-1 1.13e-5 2.27e-5 8.90e-2 7.75e-6 8.45e-6 4.43e-2I 4.01e-1 6.88e-6 1.16e-5 5.99e-3 2.66e-6 3.77e-6 2.64e-3F 3.16e-1 4.93e-6 8.90e-6 1.56e-3 1.41e-6 2.35e-6 7.35e-4RT0-P0 C 8.69e-1 1.18e-5 2.12e-5 1.06e-1 6.60e-6 5.87e-6 5.40e-2I 4.33e-1 6.91e-6 1.16e-5 5.83e-3 2.62e-6 3.69e-6 2.56e-3F 3.35e-1 5.00e-6 8.83e-6 1.54e-3 1.40e-6 2.32e-6 7.30e-4MVEM C 6.54e-1 1.16e-5 2.23e-5 8.25e-2 1.07e-5 1.29e-5 4.04e-2I 4.32e-1 6.83e-6 1.15e-5 6.15e-3 2.84e-6 3.88e-6 2.70e-3F 3.35e-1 4.92e-6 8.94e-6 1.59e-3 1.47e-6 2.39e-6 7.55e-4 † C = Coarse, I = Intermediate, and F = Fine. impermeable barriers [38], rather than a flaw related to the estimates.On a more qualitative note, Figure 5 suggests that the larger matrix errors areconcentrated around the fracture tips and fracture intersections. This comes as nosurprise, since these are exactly where singularities are located [14].
For our last numerical application, we chose benchmark 2.1 from [10]. The domainconsists of nine intersecting fractures embedded in a unit cube, as shown in the leftFigure 6. This results in an intricate network with 106 subdomains and 270 interfacesof different dimensionality.The flow inside the domain is driven by inlet (pink upper corner) and outlet (pur-ple lower corner) boundary conditions as shown in the left Figure 6. The benchmarkassigns heterogeneous permeability to the matrix subdomain, whereas the fracturesare assumed to be highly conductive. For the exact geometric setup, material con-stants, and the complete description of the problem, we refer to [10].As we did in section 6.1, we collect the local errors of subdomains and interfacesof equal dimensionality. We summarize the results in Table 4. As in the previouscases, we have local and global convergence for all four numerical methods. Again,MPFA, RT0-P0, and MVEM show very similar results, while TPFA show larger22rrors. Finally, in the right Figure 6, we show the error estimates for the wholefracture network, where it becomes evident that the errors are concentrated close tothe inlet boundary.
In this paper, we derived a posteriori error estimates for mixed-dimensional ellipticequations, which, in particular, model the incompressible single-phase flow in frac-tured porous media. Based on a functional approach, we derived abstract and fullycomputable upper bounds with no undetermined constants when locally-mass numer-ical methods are employed.To the best of our knowledge, our estimators are the first of their kind to providean upper bound for the error (of the primal and dual variables) in the case of ageneric fracture network without severe limitations in terms of the coupling betweenthe grids.Finally, our estimates have been thoroughly tested with numerical approximationsobtained with classical locally-mass conservative methods, namely TPFA, MPFA,RT0-P0, and MVEM. In particular, we perform two- and three-dimensional valida-tions with non-matching and matching grids. Moreover, we applied our estimatorsto two- and three-dimensional community benchmark problems exhibiting challeng-ing fracture networks and large material contrasts. In both cases, we have obtainedsatisfactory results.
References [1] E. Ahmed, S. Ali Hassan, C. Japhet, M. Kern, and M. Vohral´ık,
A posteriorierror estimates and stopping criteria for space-time domain decomposition fortwo-phase flow between different rock types , SMAI J. Comput. Math. (2019),195–227. MR 4062386[2] E. Ahmed, J. Jaffr´e, and J. E. Roberts, A reduced fracture model for two-phaseflow with different rock types , Math. Comput. Simulation (2017), 49–70. MR3624734[3] E. Ahmed, J. M. Nordbotten, and F. A. Radu,
Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iter-ative schemes for coupled poromechanics problems , J. Comput. Appl. Math. (2020), 112312, 25. MR 3983155[4] E. Ahmed, F. A. Radu, and J. M. Nordbotten,
Adaptive poromechanics com-putations based on a posteriori error estimates for fully mixed formulations ofBiot’s consolidation model , Comput. Methods Appl. Mech. Engrg. (2019),264–294. MR 3899054[5] C. Alboin, J. Jaffr´e, J. E. Roberts, and C. Serres,
Modeling fractures as interfacesfor flow and transport in porous media , Fluid flow and transport in porous media:mathematical and numerical treatment (South Hadley, MA, 2001), Contemp.Math., vol. 295, Amer. Math. Soc., Providence, RI, 2002, pp. 13–24. MR 1911534[6] S. S. Antman,
Nonlinear problems of elasticity , Applied Mathematical Sciences,vol. 107, Springer-Verlag, New York, 1995. MR 1323857237] I. Babuˇska and W. C. Rheinboldt,
Error estimates for adaptive finite elementcomputations , SIAM J. Numer. Anal. (1978), no. 4, 736–754. MR 483395[8] M. Bebendorf, A note on the Poincar´e inequality for convex domains , Z. Anal.Anwendungen (2003), no. 4, 751–756. MR 2036927[9] Z. Belhachmi, A posteriori error estimates for the 3D stabilized mortar finiteelement method applied to the Laplace equation , M2AN Math. Model. Numer.Anal. (2003), no. 6, 991–1011. MR 2026405[10] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gl¨aser, E. Keilegavlen,A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, S. Burbulla, P. Devloo, O. Du-ran, M. Favino, J. Hennicker, I-H. Lee, K. Lipnikov, R. Masson, K. Mosthaf,M. G. C. Nestola, C-F. Ni, K. Nikitin, P. Sch¨adle, D. Svyatskiy, R. Yanbarisov,and P. Zulian, Verification benchmarks for single-phase flow in three-dimensionalfractured porous media , Advances in Water Resources (2021), 103759.[11] D. Boffi, F. Brezzi, and M. Fortin,
Mixed finite element methods and applications ,Springer Series in Computational Mathematics, vol. 44, Springer, Heidelberg,2013. MR 3097958[12] W. M. Boon and J. M. Nordbotten,
Stable mixed finite elements for linear elas-ticity with thin inclusions , Computational Geosciences (2020), 1–18.[13] W. M. Boon, J. M Nordbotten, and J. E. Vatne,
Functional analysis and exte-rior calculus on mixed-dimensional geometries , Annali di Matematica Pura edApplicata (1923-) (2020), 1–33.[14] W. M. Boon, J. M. Nordbotten, and I. Yotov,
Robust discretization of flow infractured porous media , SIAM J. Numer. Anal. (2018), no. 4, 2203–2233. MR3829517[15] F. Brezzi, Richard S. Falk, and L. Donatella Marini, Basic principles of mixedvirtual element methods , ESAIM Math. Model. Numer. Anal. (2014), no. 4,1227–1240. MR 3264352[16] H. Chen and S. Sun, A residual-based a posteriori error estimator for single-phaseDarcy flow in fractured porous media , Numer. Math. (2017), no. 3, 805–839.MR 3660303[17] P. G. Ciarlet,
Mathematical elasticity. Vol. II , Studies in Mathematics and itsApplications, vol. 27, North-Holland Publishing Co., Amsterdam, 1997, Theoryof plates. MR 1477663[18] S. Cochez-Dhondt, S. Nicaise, and S. I. Repin,
A posteriori error estimates forfinite volume approximations , Math. Model. Nat. Phenom. (2009), no. 1, 106–122. MR 2483555[19] L. B da Veiga, F. Brezzi, L. D. Marini, and A. Russo, Mixed virtual elementmethods for general second order elliptic problems on polygonal meshes , ESAIMMath. Model. Numer. Anal. (2016), no. 3, 727–747. MR 3507271[20] C. D’Angelo and A. Quarteroni, On the coupling of 1D and 3D diffusion-reactionequations. Application to tissue perfusion problems , Math. Models Methods Appl.Sci. (2008), no. 8, 1481–1504. MR 24398472421] C. D’Angelo and A. Scotti, A mixed finite element method for darcy flow in frac-tured porous media with non-matching grids , ESAIM: Mathematical Modellingand Numerical Analysis - Mod´elisation Math´ematique et Analyse Num´erique (2012), no. 2, 465–489 (en).[22] A. Ern and M. Vohral´ık, A posteriori error estimation based on potential andflux reconstruction for the heat equation , SIAM J. Numer. Anal. (2010), no. 1,198–223. MR 2608366[23] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Ste-fansson, and A. Tatomir, Benchmarks for single-phase flow in fractured porousmedia , Advances in Water Resources (2018), 239–258.[24] L. Formaggia, A. Fumagalli, A. Scotti, and P. Ruffo,
A reduced model for Darcy’sproblem in networks of fractures , ESAIM Math. Model. Numer. Anal. (2014),no. 4, 1089–1116. MR 3264347[25] A. Fumagalli and E. Keilegavlen, Dual virtual element methods for discrete frac-ture matrix models , Oil & Gas Science and Technology–Revue d’IFP Energiesnouvelles (2019), 41.[26] F. Hecht, Z. Mghazli, I. Naji, and J. E. Roberts, A residual a posteriori errorestimators for a model for flow in porous media with fractures , J. Sci. Comput. (2019), no. 2, 935–968. MR 3968997[27] E. Hodneland, E. Hanson, O. Sævareid, G. Nævdal, A. Lundervold,V. ˇSolt´eszov´a, A. Z. Munthe-Kaas, A. Deistung, J. R. Reichenbach, and J. M.Nordbotten, A new framework for assessing subject-specific whole brain circula-tion and perfusion using mri-based measurements and a multi-scale continuousflow model , PLoS computational biology (2019), no. 6, e1007073.[28] E. Hodneland, X. Hu, and J. M. Nordbotten, Well-posedness, discretization andpreconditioners for a class of models for mixed-dimensional problems with highdimensional gap , arXiv preprint arXiv:2006.12273 (2020), 1–27.[29] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela,and I. Berre,
Porepy: An open-source software for simulation of multiphysicsprocesses in fractured porous media , Computational Geosciences (2021), no. 1,243–265.[30] T. Koch, K. Heck, N. Schr¨oder, H. Class, and R. Helmig, A new simulationframework for soil–root interaction, evaporation, root growth, and solute trans-port , Vadose zone journal (2018), no. 1, 1–21.[31] T. K¨oppl, E. Vidotto, and B. Wohlmuth, A local error estimate for thePoisson equation with a line source term , Numerical mathematics and ad-vanced applications—ENUMATH 2015, Lect. Notes Comput. Sci. Eng., vol. 112,Springer, [Cham], 2016, pp. 421–429. MR 3706388[32] V. Martin, J. Jaffr´e, and J. E. Roberts,
Modeling fractures and barriers as inter-faces for flow in porous media , SIAM J. Sci. Comput. (2005), no. 5, 1667–1691.MR 2142590[33] Z. Mghazli and I. Naji, Guaranteed a posteriori error estimates for a fracturedporous medium , Math. Comput. Simulation (2019), 163–179. MR 39801942534] J.-C. N´ed´elec,
Mixed finite elements in R , Numer. Math. (1980), no. 3,315–341. MR 592160[35] P. Neittaanm¨aki and S. I. Repin, Reliable methods for computer simulation , Stud-ies in Mathematics and its Applications, vol. 33, Elsevier Science B.V., Amster-dam, 2004, Error control and a posteriori estimates. MR 2095603[36] J. M. Nordbotten,
Mixed-dimensional models for real-world applications , Snap-shots of Modern Mathematics from Oberwolfach (2019), 11.[37] J. M. Nordbotten and W. M. Boon,
Modeling, structure and discretization ofhierarchical mixed-dimensional partial differential equations , Domain decompo-sition methods in science and engineering XXIV, Lect. Notes Comput. Sci. Eng.,vol. 125, Springer, Cham, 2018, pp. 87–101. MR 3989858[38] J. M. Nordbotten, W. M. Boon, A. Fumagalli, and E. Keilegavlen,
Unified ap-proach to discretization of flow in fractured porous media , Comput. Geosci. (2019), no. 2, 225–237. MR 3941939[39] J. M. Nordbotten and M. A. Celia, Geological storage of co2: modeling approachesfor large-scale simulation , John Wiley & Sons, 2011.[40] J. T. Oden and S. Prudhomme,
Goal-oriented error estimation and adaptivityfor the finite element method , Comput. Math. Appl. (2001), no. 5-6, 735–756.MR 1822600[41] L. E. Payne and H. F. Weinberger, An optimal Poincar´e inequality for convexdomains , Arch. Rational Mech. Anal. (1960), 286–292 (1960). MR 117419[42] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd orderelliptic problems , Mathematical aspects of finite element methods (Proc. Conf.,Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975), 1977, pp. 292–315. LectureNotes in Math., Vol. 606. MR 0483555[43] S. I. Repin,
A posteriori error estimation for variational problems with uniformlyconvex functionals , Math. Comp. (2000), no. 230, 481–500. MR 1681096[44] , Two-sided estimates of deviation from exact solutions of uniformly ellip-tic equations , Proceedings of the St. Petersburg Mathematical Society, Vol. IX,Amer. Math. Soc. Transl. Ser. 2, vol. 209, Amer. Math. Soc., Providence, RI,2003, pp. 143–171. MR 2018375[45] ,
A posteriori estimates for partial differential equations , Radon Series onComputational and Applied Mathematics, vol. 4, Walter de Gruyter GmbH &Co. KG, Berlin, 2008. MR 2458008[46] S. I. Repin, S. Sauter, and A. Smolianski,
Two-sided a posteriori error estimatesfor mixed formulations of elliptic problems , SIAM J. Numer. Anal. (2007),no. 3, 928–945. MR 2318795[47] T. H. Sandve, I. Berre, and J. M. Nordbotten, An efficient multi-point flux ap-proximation method for discrete fracture–matrix simulations , Journal of Compu-tational Physics (2012), no. 9, 3784 – 3800.[48] J. Varela,
Package mdestimates v.1.0, doi: 10.5281/zenodo.4453350 , January2021. 2649] R. Verf¨urth,
A review of a posteriori error estimation techniques for elasticityproblems , Comput. Methods Appl. Mech. Engrg. (1999), no. 1-4, 419–440,New advances in computational methods (Cachan, 1997). MR 1665355[50] M. Vohral´ık,
A posteriori error estimates for finite volume and mixed finite ele-ment discretizations of convection-diffusion-reaction equations , Paris-Sud Work-ing Group on Modelling and Scientific Computing 2006–2007, ESAIM Proc.,vol. 18, EDP Sci., Les Ulis, 2007, pp. 57–69. MR 2404896[51] ,
Unified primal formulation-based a priori and a posteriori error analysisof mixed finite element methods , Math. Comp. (2010), no. 272, 2001–2032.MR 2684353[52] M. F. Wheeler and I. Yotov, A posteriori error estimates for the mortar mixedfinite element method , SIAM J. Numer. Anal. (2005), no. 3, 1021–1042. MR2177794[53] B. I. Wohlmuth, Hierarchical a posteriori error estimators for mortar finite ele-ment methods with Lagrange multipliers , SIAM J. Numer. Anal. (1999), no. 5,1636–1658. MR 1710465[54] O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedurefor practical engineering analysis , Internat. J. Numer. Methods Engrg. (1987),no. 2, 337–357. MR 875306 27 Appendix
For completeness, we will write the exact expressions for the pressure, velocities, andsource terms for the numerical validations presented in section 5.We will conveniently define the following quantities for notational compactness: α ( x ) = x − ,β ( x ) = x − , β ( x ) = x − ,γ ( x ) = x − , γ ( x ) = x − . Exact forms for validation from section 5.1
The exact pressure in the bulk p is given by the shortest distance from any pointin Ω to the fracture Ω . Hence, p can be written as a piecewise function withthree different regions, namely bottom, middle, and top. We can therefore write thepressures for the bulk and fracture as follow p = (cid:0) α + β (cid:1) / , < x < , (cid:0) α (cid:1) / , ≤ x ≤ , (cid:0) α + β (cid:1) / , < x < , (8.1a) p = − . (8.1b)The tangential velocities are obtained via Darcy’s law, which results in u = (cid:20) − α ( α + β ) / , − β ( α + β ) / (cid:21) T , < x < (cid:20) − ( α ) / α , (cid:21) T , ≤ x ≤ , (cid:20) − α ( α + β ) / , − β ( α + β ) / (cid:21) T , < x < , (8.2a) u = 0 . (8.2b)The exact source terms are obtained using the mass conservation equation: f = α + β ( α + β ) / − ( α + β ) / , < x < , , ≤ x ≤ , α + β ( α + β ) / − ( α + β ) / , < x < , (8.3a) f = − . (8.3b)28 xact forms for validation from section 5.2 Analogously to the two-dimensional case, the pressure p inside Ω is given by theshortest distance from any point in Ω to the fracture Ω . This gives rise to a piecewisefunction with nine different regions. The exact forms for the pressure inside the bulkand fracture are given by p = (cid:0) α + β + γ (cid:1) / , < x < , < x < , (cid:0) α + β (cid:1) / , < x < , ≤ x ≤ , (cid:0) α + β + γ (cid:1) / , < x < , < x < , (cid:0) α + γ (cid:1) / , ≤ x ≤ , < x < , (cid:0) α (cid:1) / , ≤ x ≤ , ≤ x ≤ , (cid:0) α + γ (cid:1) / , ≤ x ≤ , < x < , (cid:0) α + β + γ (cid:1) / , < x < , < x < , (cid:0) α + β (cid:1) / , < x < , ≤ x ≤ , (cid:0) α + β + γ (cid:1) / , < x < , < x < , (8.4a) p = − . (8.4b)The exact tangential velocities are given by u = (cid:20) − α ( α + β + γ ) / , − β ( α + β + γ ) / , − γ ( α + β + γ ) / (cid:21) T , < x < , < x < , (cid:20) − α ( α + β ) / , − β ( α + β ) / , (cid:21) T , < x < , ≤ x ≤ , (cid:20) − α ( α + β + γ ) / , − β ( α + β + γ ) / , − γ ( α + β + γ ) / (cid:21) T , < x < , < x < , (cid:20) − α ( α + γ ) / , , − γ ( α + γ ) / (cid:21) T , ≤ x ≤ , < x < , (cid:20) − ( α ) / α , , (cid:21) T , ≤ x ≤ , ≤ x ≤ , (cid:20) − α ( α + γ ) / , , − γ ( α + γ ) / (cid:21) T , ≤ x ≤ , < x < , (cid:20) − α ( α + β + γ ) / , − β ( α + β + γ ) / , − γ ( α + β + γ ) / (cid:21) T , < x < , < x < , (cid:20) − α ( α + β ) / , − β ( α + β ) / , (cid:21) T , < x < , ≤ x ≤ , (cid:20) − α ( α + β + γ ) / , − β ( α + β + γ ) / , − γ ( α + β + γ ) / (cid:21) T , < x < , < x < , (8.5a) u = [0 , T . (8.5b)29inally, the source terms are given by f = α + β + γ ( α + β + γ ) / − ( α + β + γ ) / , < x < , < x < , α + β ( α + β ) / − ( α + β ) / , < x < , ≤ x ≤ , α + β + γ ( α + β + γ ) / − ( α + β + γ ) / , < x < , < x < , α + γ ( α + γ ) / − ( α + γ ) / , ≤ x ≤ , < x < , , ≤ x ≤ , ≤ x ≤ , α + γ ( α + γ ) / − ( α + γ ) / , ≤ x ≤ , < x < , α + β + γ ( α + β + γ ) / − ( α + β + γ ) / , < x < , < x < , α + β ( α + β ) / − ( α + β ) / , < x < , ≤ x ≤ , α + β + γ ( α + β + γ ) / − ( α + β + γ ) / , < x < , < x < , (8.6a) f = − ..