Abstract
Mixed-dimensional elliptic equations exhibiting a hierarchical structure are commonly used to model problems with high aspect ratio inclusions, such as flow in fractured porous media. We derive general abstract estimates based on the theory of functional a posteriori error estimates, for which guaranteed upper bounds for the primal and dual variables and two-sided bounds for the primal–dual pair are obtained. We improve on the abstract results obtained with the functional approach by proposing four different ways of estimating the residual errors based on the extent the approximate solution has conservation properties, i.e., (1) no conservation, (2) subdomain conservation, (3) grid-level conservation, and (4) exact conservation. This treatment results in sharper and fully computable estimates when mass is conserved either at the grid level or exactly, with a comparable structure to those obtained from grid-based a posteriori techniques. We demonstrate the practical effectiveness of our theoretical results through numerical experiments using four different discretization methods for synthetic problems and applications based on benchmarks of flow in fractured porous media.
1 Introduction
Mixed-dimensional partial differential equations (mD-PDEs) arise when partial differential equations interact on domains of different topological dimensions [43]. Prototypical examples include models of thin inclusions in elastic materials [8, 17, 22], blood flow in human vasculature [25, 33, 37], root water uptake systems [36], and flow in fractured porous media [3, 7, 30]. The latter example has an appealing mathematical structure, in that the model equations allow for a hierarchical representation where each subdomain (matrix, fractures, fracture intersections, and intersection points) only has direct interaction with subdomains of topological dimension one higher or one lower [19]. Such hierarchical mD-PDEs are the topic of the current paper.
![Fig. 1
Example geometries falling within the context of hierarchical mixed-dimensional geometries studied herein. Left panel corresponds to a 2D benchmark problem [29] while the two remaining correspond to 3D benchmark problems [14].](/document/doi/10.1515/jnma-2022-0038/asset/graphic/j_jnma-2022-0038_fig_001.jpg)
mD-PDEs are intrinsically linked to the underlying geometric representation, which, in a certain sense, generalizes the usual notion of the domain. One can then define sets of suitable functions (and function spaces) on this geometry, and these sets are then naturally interpreted as mixed-dimensional (mD) functions. Exploiting this concept, one can generalize the standard differential operators to mappings between mD functions and thus obtain an mD calculus. The fact that this mD calculus inherits standard properties of calculus, particularly partial integration (relative to suitable inner products), a de Rham complex structure, and a Poincaré–Friedrichs inequality, was recently established using the language of exterior calculus on differential forms [18].
The inherent geometric generality of hierarchical mD-PDEs also demand the same level of abstraction of a posteriori error estimation techniques. This requirement makes error estimates of the functional type particularly well-suited for the task [42, 47, 55, 56, 57, 59]. The most attractive feature of this approach is that error estimates are derived using purely functional methods [57]. The bounds are therefore agnostic to the way approximated solutions are obtained in the energy space, and the only undetermined constants arise from Poincaré-type inequalities [38].
However, unlike other types of error estimates [6, 46, 61, 63, 66], this generality makes standard functional estimates of limited applicability to hierarchical elliptic mD-PDEs due to the following reasons: (1) for general fracture networks, the mixed-dimensional Poincaré constant is not easily computable, and (2) since Poincaré constants are proportional to the diameter of the physical domain, residual estimators cannot exhibit superconvergent properties.
To circumvent the aforementioned issues, we exploit the fact that Poincaré-type inequalities imply weighted norms [50, 53], and use spatially-dependent weights to control the residual norms. We show both theoretically and numerically that this treatment leads to sharper estimates when approximations to the exact solution satisfy mass conservation in a given partition of the domain.
In view of the preceding discussion, our aim is therefore to obtain a posteriori error estimates for the approximate solution to the mD scalar elliptic equation [18, 19, 44], where the mD Laplace equation for geometries such as those illustrated in Fig. 1, is described in detail in Section 3.
We remark that while a broad range of a posteriori error techniques are available for mono-dimensional problems, existing error bounds for mD models are far more scarce. Moreover, the ones available, are restricted to specific cases (e.g., in the context of mortar methods [13, 52, 64, 65] and fractured porous media [21, 32, 40]) with far less geometric generality than what we present here. Thus, for practical problems, a posteriori error bounds for mD geometries have until now essentially not been available.
The rest of the paper is structured as follows: Section 2 is devoted to introducing the model problem, functional spaces, and variational formulations for the case of a single 1D fracture embedded in a 2D matrix. The section is concluded by providing a first upper bound for the primal variable. In Section 3, we generalize the results from Section 2 to the case of fracture networks and introduce the necessary tools to perform the a posteriori analysis in an mD setting. After reviewing necessary tools from functional analysis in Section 4, in Section 5, we provide our main results starting from a generic abstract estimate and then considering specific cases depending upon the degree of accuracy at which residual terms are approximated. In Section 6, we introduce the approximated problem using mixed-finite element methods and thus make the estimates concrete. Section 7 deals, respectively, with numerical validations and practical applications of the derived bounds. Finally, in Section 8, we present our concluding remarks.
2 Upper bounds for a single fracture
In this section, we introduce the model problem together with functional spaces and the variational formulations for the case of a single 1D line embedded in a 2D matrix, as illustrated in Fig. 2. Furthermore, a first upper bound for the primal variable is derived following the classical functional approach. We remark that the case of a single fracture embedded in a matrix has been analyzed before. For example, [21] and [32] proposed error estimators based on the residual approach, whereas [40] obtained guaranteed a posteriori error estimates using the approach of Vohralík [63].
2.1 The model problem for a single fracture
Before writing the set of equations describing general fracture networks, let us first introduce the governing equations of a simpler configuration; that is, a unit square domain Y ⊂ ℝ2 decomposed as a 1D fracture Ω1 embedded in a 2D matrix Ω2, as shown in the left panel of Fig. 2. Interfaces Γ1 and Γ2, at each side of Ω1, establish the link between Ω2 and Ω1. The model presented below is well-established for these problems, and we point the reader to the references for further justification of this system [19, 39, 44].
A horizontal 1D fracture embedded in a 2D matrix. Left: Subdomains and interfaces. Right: Boundary conditions. For the fracture, the purple square denotes a no-flux boundary condition, whereas the green square a Dirichlet boundary condition. Note that ∂1Ω2, Γ1, Ω1, Γ2, ∂2Ω2, all coincide spatially. For illustrative purposes, however, they are placed in different locations.
The strong form of the governing equations in Ω2 reads
Here, (2.1a) is the mass conservation equation, u2 is the matrix velocity, and f2 an external source. The fluid velocity is given by the standard Darcy’s law (2.1b), where 𝓚2 is the matrix permeability; a bounded, symmetric, and positive-definite 2 × 2 tensor, and p2 is the fluid pressure.
Equations (2.1c) and (2.1d) require that at each side of the internal boundary of Ω2, the normal component of u2 to match the interface (mortar) fluxes λ1 and λ2. To fix the direction of the normal vector on internal boundaries, we require n2 pointing from the higher- to the lower-dimensional subdomain. No flux conditions are prescribed in (2.1e), where u2⋅n2 represents the outer normal flux across ∂NΩ2. Finally, Dirichlet boundary conditions are imposed in (2.1f), where gD,2 is a prescribed function on the Dirichlet boundary.
In the fracture Ω1, the equations are given by
In (2.2a),
The tangential velocity u1 is again expressed via Darcy’s law (2.2b), where in a slight abuse of notation, we use 𝓚1 to refer to the tangential component of the fracture permeability, which is again assumed to be positive and bounded from above. Finally, (2.2c) and (2.2d) are the Neumann and Dirichlet boundary conditions, respectively. Again, we use gD,1 to denote a prescribed function on the Dirichlet part of the fracture boundary.
To close the system of equations, we must specify a constitutive relationship for the interface fluxes. Here, we use a Darcy-type law [39], where mortar fluxes are linearly related to pressure jumps
with κ1 and κ2 representing the effective normal permeability on Γ1 and Γ2, respectively. We restrict our analysis to the case where κ1 and κ2 are non-degenerate. Thus, following [19], we further require the existence of two constants γ1 and γ2 such that 0 < γ1 ⩽
2.2 Functional spaces and variational formulations
Let us now present the primal weak formulation of the single fracture model from the previous section. To this aim, consider first the energy space with vanishing traces on Dirichlet boundaries
and the product spaces
Furthermore, let 〈⋅, ⋅〉Ωi and 〈⋅, ⋅〉Γj denote respectively the L2–inner products on Ωi and Γj, and ∥⋅∥Ωi and ∥⋅∥Γj the relevant L2–norms. Finally, we denote by g = [g1, g2] ∈ H1(Ω) two functions extending the boundary data into the domains, and thus satisfying tr∂DΩi gi = gD,i. We now state the primal weak problem as the following.
Definition 2.1
(primal weak formulation for a single fracture). Let p = [p1, p2] and g = [g1, g2] ∈ H1(Ω). Then find p ∈
Refer to Appendix A.1 for the derivation of the primal weak form from the strong form in Section 2.1. We see directly from equation (2.6) that the primal weak form has a minimization structure subject to the stated conditions on 𝓚i and κj, and well-posedness follows by standard arguments.
A dual weak form for the model problem, with explicit representation of the subdomain velocities and mortar fluxes, can also be constructed. We first define the space H(div; Ωi, ∂X Ω) as the space of L2-vector functions on Ωi with weak divergence in L2(Ωi) and zero trace on the part of the boundary indicated by ∂X Ω. Then, we denote the product spaces of H(div)-functions that are zero on Neumann, and on Neumann and internal boundaries as:
Furthermore, we define the L2-product spaces on the domains:
With these spaces in hand, we consider the standard linear extension operator from internal boundaries onto domains denoted 𝓡j := L2(Γj) → H(div; Ω2, ∂N Ω2), such that 𝓡j satisfies for all λj ∈ L2(Γj)
The precise choice of the extension operator 𝓡j is not important; however, the natural choice based on the solution of an auxiliary elliptic equation is reasonable [19]. We naturally extend the definition of 𝓡j to 𝓡 := L2(Γ) → V by requiring that for [λ1, λ2] ∈ L2(Γ), then [u1, u2] = 𝓡λ satisfies u1 = 0 and u2 = 𝓡1 λ1 + 𝓡2 λ2.
The above constructions allow us to represent subdomain fluxes as
where u0 ∈ V0 and λ ∈ L2(Γ). This motivates the construction of a compound H(div)-type spaces, as
This construction will become key when we generalize to more complex geometries in the next section.
Remark 2.1
(on the regularity of H(div; Ω, Γ)). It is worth remarking that the restriction of the space H(div; Ω, Γ) to the domain Ω2 has slightly enhanced regularity relative to the standard space H(div; Ω2), as this latter space has normal traces which do not lie in L2(Γ1) nor L2(Γ2).
Definition 2.2
(dual weak formulation for a single fracture). Let u0 = [u0, 1, u0, 2], λ = [λ1, λ2], p = [p1, p2]. Then find (u0, λ, p) ∈ H(div;Ω, Γ) × L2(Ω) such that
Refer to Appendix A.2 for the derivation.
Remark 2.2
(well-posedness). The variational formulation from Definition 2.2 can be classified as a saddle point structure, for which well-posedness results have been established for fracture networks (see, e.g., [19, Th. 2.5]).
2.3 A first a posteriori error estimate for the primal variable
Having the functional spaces and weak formulations formally introduced, in this section, we provide a first upper bound for an approximation to the primal variable q = [q1, q2] ∈
Theorem 2.1
(a first upper bound for the primal variable). Let p ∈
with
where CΩ1 and CΩ2 are the permeability-weighted Poincaré–Friedrichs constants for Ω1 and Ω2:
Proof
Refer to Appendix B for the proof.□
Remark 2.3
(nature of the estimators). The upper bound (2.15) is a guaranteed upper bound for the deviation between the primal solution p ∈
Remark 2.4
(sharpness of the estimates). The estimates above are in principle sharp, as can be shown by standard arguments [57]. However, in practice, we will often have access to additional information about the approximate solution (most commonly if it is derived with a local conservation property). This allows for improvements in the residual estimators (2.16f) and (2.16e), as we will show in Section 5.2.
It is clear that even for this fairly simple configuration, the variational formulations (and the analysis in general) can be quite cumbersome. The situation escalates in complexity when intersecting fractures (see Fig. 3) are part of the geometric configuration, in particular as the proof of Theorem 2.1 relies on all subdomains having some non-vanishing Dirichlet boundary. Indeed, when floating subdomains (e.g., fully embedded fractures or isolated rock domains) are present in the fracture network, the standard procedure used in Theorem 2.1 can no longer be applied directly. Thus, in the remainder of the paper, we deal with these challenges in a more general framework.
Mixed-dimensional geometric decomposition of a fracture network. Left: The domain Y is decomposed into two 2D matrices (Ω9 and Ω10), four 1D fractures (Ω5, Ω6, Ω7, and Ω8), one 0D fracture intersection point (Ω4), and three 0D fracture end-points (Ω1, Ω2, Ω3). Note that we allow fractures and other lower-dimensional subdomains to form parts of the boundary of the domain (e.g., Ω5 with its endpoints Ω1 and Ω2). Center: Interfaces between subdomains. Right: Subdomain boundaries. Internal boundaries are depicted in red, whereas fracture’s boundaries touching the ambient boundary are depicted in green.
3 Extension to fracture networks
In this section, we extend the single fracture model to account for several subdomains as part of a general fracture network. Our vocabulary is motivated by the physical case of n = 3, where the surrounding rock is composed of simply connected 3D subdomains, fractures are simply connected planar 2D subdomains, the intersection between such fractures are 1D lines, and the intersection between fracture intersections are 0D points (see Fig. 3 for an example with n = 2).
We start with the classical description and then introduce the mD notation. The rest of the section is devoted to introducing key tools that are necessary to perform the analysis in an mD setting.
3.1 Mixed-dimensional geometric representation
The derivation of a posteriori estimates for generic fracture networks greatly benefits from an mD decomposition of the domain of interest, and we therefore follow the approach of [19]. We start by considering an n–dimensional contractible domain Y ⊂ ℝn, n ∈ {2, 3}, decomposed into m planar, open and non-intersecting subdomains Ωi of different dimensionality di = d(i), such that Y =
We adopt a structure where neighboring subdomains one dimension apart are connected via interfaces, denoted by Γj for j ∈ {1, …, M}. To be precise, let Γj be the interface between subdomains indexed by ǰ and ĵ of dimension d and d + 1, respectively. Then Γj = Ωǰ (see center panel of Fig. 3), and furthermore, we denote the adjacent boundary of Ωĵ by Γj = ∂j Ωĵ. We emphasize that while the internal boundary ∂jΩĵ is defined to spatially coincide with the interface Γj, which in turn coincides with the lower-dimensional subdomain Ωǰ, their distinction is crucial to define variables properly.
To keep track of the connections from subdomains to interfaces, we introduce the sets 𝓢̂i and 𝓢̌i, containing the indices of the higher-dimensional (respectively lower-dimensional) neighboring interfaces of Ωi, as illustrated in the right panel of Fig. 3. These sets are dual to ǰ and ĵ defined in the previous paragraph, thus for all j ∈ 𝓢̂i, it holds that ǰ = i, while for all j ∈ 𝓢̌i, it holds that ĵ = i.
We will be interested in defining functions on the above stated partition of the domain and the interfaces. This motivates us to define the disjoint unions
A complete mixed-dimensional partitioning, including both subdomains and interfaces, is given by Ω ⊔ Γ.
In order to speak of boundary conditions, we introduce the decomposition of the boundary of Ω. Let ∂Ω be partitioned into its Neumann, Dirichlet, and internal parts. That is, we define ∂Ω = ∂NΩ ∪ ∂DΩ ∪ ∂IΩ, where
3.2 The model problem for a fracture network
Let us now present the model problem valid for m subdomains of dimensionality 0 to n, and M interfaces of dimensionality 0 to n − 1. Our model summarizes the derivations given in recent literature [19, 34, 44]. For all domains Ωi, we consider a scalar pressure pi together with a flux ui in the tangent space of the domain. On all interfaces Γj, we consider a scalar coupling flux λj, oriented as positive for flow from the higher dimensional domain Ωĵ. We will, in this section, assume sufficient regularity that the strong form makes sense, and return to the weak formulation in later sections. The governing equations from the previous section then generalize as
In (3.2a), the summation captures the contribution of fluxes from the adjacent interfaces to Ωi, and can be seen as a generalization of the second term in (2.2a). Note that for di = n, the set 𝓢̂i = ∅, and thus the jump operator, evaluates to zero in the highest-dimensional domains. Conversely, in (3.2a), the differential term ∇i ⋅ ui is void whenever di = 0, as there is no tangent space to a point in all subdomains, and indeed, we will not consider the ui defined on these domains, which justifies why equation (3.2b) are not applied to 0D domains.
We are now ready to recast the model problem in mD notation, building on the product space structures introduced in Section 2.2. Let us start by defining the mD pressure as the ordered collection of subdomain pressures 𝔭 := [pi] ∈ CΩ, i.e., scalar functions on Ω. We now decompose the fluxes as in (2.11), so that
such that u0,i satisfies u0,i ⋅ ni = 0 for all j ∈ 𝓢̌i, and where the reconstruction operator is generalized as 𝓡j : CΓj → CΩĵ satisfying:
This allows us to define the mD flux as the internal (tangential) domain fluxes and (normal) interface fluxes 𝔲 := [u0,i, λj] ∈ C0 TΩ × CΓ, i.e., the pairing of sections of the tangent bundle TΩ together with scalar functions on Γ. By the subscript C0 TΩ, we indicate that both ui ⋅ ni = 0 on all ∂jΩi, where j ∈ 𝓢̌i, and also ui ⋅ ni = 0 on ∂NΩi.
We now define a generalized divergence operator 𝔇⋅(⋅): C0 TΩ × CΓ → CΩ which acts on the mD flux in accordance with the left-hand side of (3.2a):
where 𝔮 = [qi] ∈ CΩ is a scalar function for each domain Ωi, defined by
Similarly, we define an mD gradient operator 𝔻 (⋅): CΩ → CTΩ × CΓ acting on the mD pressure in accordance with the right-hand sides of equations (3.2b) and (3.2c):
where 𝔳 = [v0,i, νj] ∈ CTΩ × CΓ has the same structure as the mD flux (but without the boundary conditions), such that for all i ∈ {1, …, m} and j ∈ {1, …, M}, it holds that
Recalling that the full flux vi is recovered from equation (3.3), we note that the second term above is simply the gradient on each subdomain. We will, in Section 3.3, further justify the terminology ‘divergence’ and ‘gradient’ due to the fact that these operators satisfy an integration-by-parts property with respect to the suitable inner products, and are thus adjoints (subject to appropriate boundary conditions).
Material parameters are collected into the mD permeability 𝔎: CTΩ × CΓ → CTΩ × CΓ, defined such that for
then from the model given in equation (3.2), we recognize the desired relationships
The second term, corresponding to Darcy’s law, can be rewritten in terms of the decomposition 𝔲 = [u0,i, λj] from equation (3.3) as
The presence of the extra terms arising from the decomposition is analogous to that in (3.2).
We note that the restriction 𝔲 ∈ C0 TΩ × CΓ, implicitly places constraints (depending on the material constants 𝔎 and via the definition of 𝔻) on the admissible pressures 𝔭. This space of admissible pressures can be understood as the domain of the restricted operator 𝔎𝔻 : CΩ → C0 TΩ × CΓ.
In view of the mD variables and operators defined above, and subject to the right-hand side data 𝔣 = [fi] ∈ CΩ and the boundary data 𝔊D = [gD,i] ∈ C∂DΩ, a straightforward substitution of definitions shows that problem (3.2) is equivalent to the concisely stated mD elliptic problem
defined for 𝔲 ∈ C0 TΩ × CΓ and 𝔭 ∈ CΩ.
Remark 3.1
(internal Neumann boundaries). For simplicity of exposition, the domain Y is taken as contractible, and Ωi is considered a partitioning of Y. However, the reader will appreciate that these assumptions can be relaxed. Most importantly, from the perspective of applications (as discussed in Section 2.1), some internal interfaces may be modeled as impermeable, i.e., λj = 0. We refer to the remaining (permeable) interfaces as Ξ ⊂ {0, …, M}. The impermeable interfaces are then excluded from the problem, and considered as internal Neumann interfaces. To be precise, we define a reduced disjoint union of interface domains
The internal Neumann boundaries may partition the domain into disconnected parts. We refer to a subdomain as ‘Dirichlet-connected’, denoted i ∈ ξ if either (1) ∂DΩi ≠∅, or (2) there exists some j ∈ 𝓢̂i such that ĵ ∈ ξ, or (3) there exists some j ∈ 𝓢̌i such that ǰ ∈ ξ. This allows us to construct a reduced disjoint union of subdomains
All the derivations in the continuation are equally valid for these reduced product domains.
Remark 3.2
(extensions to the model equations). The results of this paper can with minor modifications be extended to non-zero Neumann boundary conditions, and with some additional effort to the class of non-planar geometries considered in [18]. However, as this generality is typically not needed for applications, we restrict the presentation as indicated above.
3.3 Variational formulations in mixed-dimensional notation
Before writing the variational formulations in mD notation, let us first define the relevant mD inner products and norms. Consider the following inner-products
and their respective induced norms
With these inner products, the previously defined mD divergence satisfy the following integration-by-parts formula [18, 19] whenever 𝔳 ∈ CTΩ × CΓ and 𝔮 ∈ CΩ:
In the above, the restriction to the boundary is denoted 𝔗X(⋅) (for X = D, N), which depending on context acts as the boundary values of pressure variables, 𝔗X(⋅): CΩ → C∂XΩ, or the normal component of flux variables, 𝔗X(⋅): CTΩ × CΓ → C∂XΩ.
From the product structure in the definition of the C and L2 spaces, the continuous spaces inherit their density from the individual subdomains to the product spaces on Ω and Γ. We can thus follow standard procedures to obtain weak extensions of the mD differential operators, the boundary restriction (trace) operators, and the corresponding function spaces [2, 10, 51]. We elaborate this below.
Due to the density of C0 TΩ × CΓ in L2 TΩ × L2Γ, the mD divergence from Section 3.2 is a densely defined unbounded linear operator on the latter space 𝔇⋅: L2Ω → L2 TΩ × L2Γ. Let us now (temporarily) use the notation (T, dom(T)) to emphasize that an operator T has domain of definition dom(T), and we denote the adjoint operator with respect to the L2 inner product by an asterisk.
We recall that the Neumann boundary is incorporated into the definition of the continuous flux spaces C0 TΩ × CΓ, thus the last term in the integration-by-parts formula (3.17), is zero. Hence, we can define a weak mD gradient and the corresponding space of weakly mD differentiable functions with zero trace on the Dirichlet boundary
Clearly, C0Ω ⊆
Considering the integration-by-parts formula again, the weak mD divergence and the corresponding space of flux functions with divergence in L2 and zero trace on the Neumann boundary H(div; Ω, Γ) can be defined as
Again C0 TΩ × CΓ ⊆ H(div; Ω, Γ), and it is appropriate to consider (𝔇⋅, H(div; Ω, Γ)) as a weak divergence. This domain of definition of the weak divergence has the interpretation of H0(div; Ωi) on all subdomains Ωi (where the subscript zero indicates zero trace on all boundaries except for Dirichlet boundaries), and L2(Γj) spaces on all interfaces Γj. Thus H(div; Ω, Γ) =
Due to the above identification of H1(Ω) and H(div; Ω, Γ) in terms of product spaces of standard function spaces on subdomains, we extend the definition of the boundary restriction operators 𝔗X(⋅) to trace operators on the weak spaces by requiring that they coincide with the standard trace operators on subdomains.
In the continuation, we will always consider the weak mD gradient and divergence, and denote these simply by 𝔻 and 𝔇⋅, respectively. Similarly, we will always consider the boundary restrictions as trace operators. The above definitions of weak mD gradient and divergence operators, and their adjoint property on the above weak spaces, has the following statements of the primal and dual weak formulations of equations (3.12) as a direct consequence.
Definition 3.1
(mixed-dimensional primal weak formulation). Let 𝔊 ∈ H1(Ω). Then find 𝔭 ∈
Definition 3.2
(mixed-dimensional dual weak formulation). Find (𝔲, 𝔭) ∈ H(div;Ω, Γ) × L2(Ω) such that
The above weak forms of the mixed-dimensional elliptic problem are well-posed for bounded coefficients [18], in the sense that there exist positive constants 𝔎0 and 𝔎∞ such that
The solutions of the primal and dual weak formulations are equivalent, and define true solutions 𝔭 ∈
4 Functional analysis tools
In this section, we summarize the main functional analysis tools we will exploit for the a posteriori analysis.
4.1 Poincaré-type inequalities
We recall the following weighted Poincaré inequalities.
Lemma 4.1
(permeability-weighted Poincaré–Friedrichs inequalities). There exist constants CΩ ⩾ CΩi ⩾ CK such that
Here, we denote by q̃Ωi and q̃K the mean value of q over the subdomain Ωi and an arbitrary di-simplex K ⊂ Ωi, respectively.
We refer to CΩ,Γ as the mixed-dimensional permeability-weighted Poincaré–Friedrichs constant (whose existence was shown in [18]), CΩi is the standard subdomain permeability-weighted Poincaré–Friedrichs constant, and CK is a local permeability-weighted Poincaré–Friedrichs constant.
It is important to mention that concrete values of CΩi are available only for a limited set of geometries (see, e.g., [20, 58, 62]). An upper bound exists for convex domains, and thus for a simplex K ⊂ Ωi we have [12, 49]:
where cK is the lower bound on the permeability within K:
The importance of this is understood if K is an element of a simplicial partition of Ωi, in which case CK scales with the mesh size hK = diam(K). This allows for superconvergent properties of residual estimators for some locally mass-conservative approximations [26, 28, 63]. We analyze these cases with further details in Section 5.2 and Remark 6.6.
4.2 Conforming flux spaces
It is often possible to verify that an approximate solution 𝔳 ∈ H(div; Ω, Γ) satisfies some stronger conservation property, that is to say, that there is some space U ⊆ L2 such that
This allows for the construction of stronger a posteriori estimates, and as such, we formalize this concept as a generalization of H(div; Ω, Γ) to ‘U-conforming flux spaces’.
Definition 4.1
(Conforming mD flux space). Let H(div; Ω, Γ; U) ⊂ H(div; Ω, Γ) be a U-conforming flux space, in the sense of
To exploit the conforming flux spaces, we must construct certain projected H1(Ω) spaces. Consider therefore U as some subspace of L2(Ω) and define U⊥ to be its orthogonal complement:
Moreover, let πU⊥ be the L2–projection onto U⊥, such that for any 𝐯 ∈ L2(Ω), πU⊥ 𝐯 ∈ U⊥ satisfies the orthogonality property:
Consider now the projected
which are defined within the class 𝓒W with unit Poincaré constants:
Indeed, such classes exist in the literature of Poincaré inequalities for weighted norms (see, e.g., [50, 53]). Note that a trivial member of 𝓒W is the inverse of the permeability-weighted mD Poincaré–Friedrichs constant μ(x) =
Remark 4.1
(on the space H(div; Ω, Γ; U)). The conforming mD flux spaces allow us to obtain sharper estimates in Section 5. However, it is important to note that the standard case U = L2(Ω) is included in our definition, for which the orthogonal complement is void, and the projection πW = I; thus W =
4.3 Bilinear forms and energy norms
For the a posteriori analysis, we will need the next two mD bilinear forms and their induced energy norms
which are related via
We also define the full norm for a mixed-dimensional pair of primal and dual variables as
Note that the last norm will depend on the eventual choice of μ−1, which we emphasize must be from the class μ ∈ 𝓒W, as defined in the preceding section.
5 A posteriori error estimates
This section is devoted to obtaining the error bounds for our model problem. First, we provide general abstract estimates, and later we focus on the evaluation of the different bounds.
5.1 General abstract estimates
Let us now present the general abstract bounds. We formalize the main results presented in Section 3 and extend the ones presented in Theorem 2.1 in the following theorem.
Theorem 5.1
(general abstract a posteriori error bounds). Let the error majorant be defined as
where
valid for all 𝔮 ∈
Let 𝔭 ∈
where
Let 𝔲 ∈ H(div; Ω, Γ) be the solution to (3.21) and 𝔳 ∈ H(div; Ω, Γ; U) be arbitrary. Then
where
Let 𝔭 ∈
where
Proof
Due to the construction of mixed-dimensional product spaces and the adjoint property of the weak differential operators, the proof from the mono-dimensional case can (to a large extent) be applied directly [59]. A notable deviation from the standard proofs is the use of conforming flux spaces, and the inclusion of the Poincaré-constants in the weights 𝓒W. The full proof is included for completeness in Appendix C.□
Remark 5.1
(non-conforming approximations). Referring again to the general setting of mD calculus, it has been shown that the differential operators form part of a cochain complex, and that an mD Helmholtz decomposition exists [18]. Thus, by realizing the above constructions as Hilbert complexes, the above error bounds can be extended also to non-conforming approximations following, e.g., [47, Th. 4.7]. However, as a main objective of our work is to obtain bounds based on conforming properties of the approximations, we will not pursue non-conforming approximations in this work.
5.2 Evaluation of the majorant
The aim of this section is to provide concrete forms of the majorant 𝓜(𝔮, 𝔳, 𝔣, μ) from Theorem 5.1 depending upon the choices of the weights μ. For this purpose, consider once again the definition of the majorant
The estimation of the first term ηDF(𝔮, 𝔳) is independent of the weights μ. Indeed, by applying (4.11), it is straightforward to see that
The terms ηDF∥,K and ηDF⊥,K measure the diffusive flux error in the tangential and normal directions associated with the subdomain element K ∈ 𝓣Ωi and the mortar element K ∈ 𝓣Γj, respectively.
To complete the evaluation of the majorant, we are left with the estimation of ηR(𝔳, 𝔣, μ), which depends on the choices of μ. Recall that this term measures the mismatch in satisfying the conservation equation in each subdomain. To be precise, there are four main types of conforming fluxes; standard L2-conforming, subdomain conservation, grid level (local) conservation, and point-wise. The quality of the residual balance can be verified explicitly before applying the a posteriori estimates, and thus is not considered an assumption in the theory. Below, we make precise the aforementioned cases.
5.2.1 No mass-conservation
Assume nothing is known about the approximation of the residual terms beyond the L2 structure. We indicate this case by the abbreviation ‘NC’, and set UNC = L2, and 𝔳 ∈ H(div; Ω, Γ; UNC) = H(div; Ω, Γ). Then
Using (4.10) and the mD Poincaré inequality (4.1a), one obtains the following bound, which is the weakest bound available within the class of bounds considered in this paper:
Here, ηR,K;NC denotes the local residual error for non-conservative approximations. The majorant when mass conservation cannot be guaranteed at any level is then given by
and it follows from the above that this is an upper bound, 𝓜 ⩽ 𝓜NC.
5.2.2 Subdomain mass-conservation
Due to the structure of the equations, where interface fluxes are stated explicitly, many approximations will have mass conservation satisfied in a subdomain level, which is in a sense a compatibility condition on the floating domains Ωi. We indicate this case by the abbreviation ‘SC’. In particular, the divergence 𝐯 = [ri] = 𝔇⋅ 𝔳 ∈ USC satisfies for all i ∈ {1, …, m} where ∂DΩi = ∅,
Thus, by definition
This case represents an improvement relative to the previous one, in the sense that we can now employ the subdomain Poincaré constants instead of the mD constant. Let us make this point precise in the following lemma.
Lemma 5.1
Let W =
Then, μ(x) =
Proof
Using the Poincaré inequality (4.1c) and the fact that the sum of broken norms is weaker than the full norm, the following result holds
□
In view of Lemma 5.1, ηR can be bounded as
where ηR,K;SC are the local residual estimators for subdomain mass-conservative approximations. The majorant for this case is given by
This estimate is sharper than that in the preceding section, since CΩi ⩽ CΩ,Γ, thus whenever the assumptions of this section are satisfied, it holds that 𝓜 ⩽ 𝓜SC ⩽ 𝓜NC.
Note that (5.13) is identical in structure to the residual estimators (2.16e) and (2.16f) obtained in Theorem 2.1. However, they are fundamentally different in the sense that (5.13) do not require all subdomains to have a non-empty Dirichlet part but rather mass to be conserved in each subdomain Ωi.
5.2.3 Local mass-conservation
By choice of numerical method, it is often easy to verify that mass is conserved on an element basis in a subdomain partition. We indicate this case by the abbreviation ‘LC’. As in the preceding section, this implies that the divergence 𝐯 = [ri] = 𝔇⋅ 𝔳 ∈ ULC then satisfies for all K ⊂ 𝓣Ωi that
where 𝓣Ωi denotes a finite partition of Ωi (typically a simplicial grid). In this case, ULC contain functions having zero mean on each element K ∈ 𝓣Ωi, and from (5.14) we see that
We will consider the slightly weaker case, where (5.14) is only required to hold for all ‘non-Dirichlet boundary’ elements, that is for all elements where ∂ K∩ ∂DΩ = ∅. This is sufficient for the results from Lemma 5.1 to be extendable to the grid level by considering the space
Lemma 5.1 now applies without modification, and weights μ(x) ⩾
where ηR,K;LC are the local residual estimators for locally mass-conservative approximations. Using the above results, the majorant for locally mass-conservative approximations reads
The local residual estimates ηR,Ω;LC correspond to the ones previously obtained by [26, 63] for mono-dimensional problems subject to a flux equilibration step. Since CK ⩽ CΩi, then as before, whenever the assumptions of this section are satisfied, it holds that 𝓜 ⩽ 𝓜LC ⩽ 𝓜SC ⩽ 𝓜NC.
Remark 5.2
(fully computable residual estimators). Unlike estimators obtained with residual methods (containing unknown constants [11, 35]) or a purely functional approach such as in Sections 5.2.1 and 5.2.2 (containing constants that are generally difficult to determine [57]), estimators such as (5.15) contain only known local constants depending on the mesh size and material parameters. This justifies the claim that these estimators are fully computable.
5.2.4 Exact mass-conservation
Methods with local mass conservation, as discussed in the previous section, when applied to problems where the RHS data 𝔣 is zero or piecewise constant, can then often be verified to have an exact (pointwise) conservation property. We indicate this case by the abbreviation ‘EC’, for which 𝔣 = 𝔇⋅𝔳, so that UEC = 0 and
This case can also be seen as the limiting case of local mass conservation for a family of grid partitions where hK → 0.
5.2.5 Summary of majorants and subdomain errors
With the obtained majorants, we can define the corresponding upper bounds for the errors of the primal, dual, and primal–dual variables.
Definition 5.1
Let α = NC, SC, LC, EC, corresponding to the flux conformity spaces Uα discussed in the preceding sections. Then, in view of the results from Theorem 5.1 and the majorants (5.9), (5.13), (5.16), and (5.17), the upper bounds for the error in the primal, dual, and primal–dual pair, for arbitrary approximations 𝔮 ∈
while the lower bound for the error in the primal–dual pair is
It is our interest not only to measure local errors, but also to distinguish between subdomain and interface errors. This motivates the definition of the following error indicators.
Definition 5.2
(subdomain and interface error indicators). Let α = EC, LC, SC, NC. Then, we will denote by εΩi;α and εΓj the subdomain and interface error indicators, defined by
We emphasize that while the majorants provide guaranteed bounds, the subdomain and interface error indicators can only be expected to correlate with the error.
6 Concrete bounds for locally mass-conservative approximations
In this section, we will make the evaluation of the bounds concrete by providing explicit approximations to (3.21) using the lowest-order mixed-finite element method (MFEM).
6.1 Grid partitions
Ultimately, a posteriori estimates are primarily applied to approximations that are defined on computational grids. We therefore, in this section, summarize the relevant notation for grids.
Let us start by defining the partitions of the domains of interest. To this aim, denote by 𝓣Ωi, 𝓣Γj, and 𝓣∂jΩi the partitions of Ωi, Γj, and ∂jΩi, respectively. Moreover, let
Here, we only consider simplicial partitions. In particular, we require all elements K ⊂ Ωi to be strictly non-overlapping simplices of dimension dK = di. We use hK to denote the diameter of K, and define hΩi = maxhK 𝓣Ωi, hΓj = maxhK 𝓣Γj, and h∂jΩi = maxhK 𝓣∂jΩi.
Left: Matching coupling between the grids 𝓣Ωĵ, 𝓣Γj, and 𝓣Ωǰ. Right: Degrees of freedom involved in the coupling between a 2D higher-dimensional cell, a 1D mortar-cell, and a 1D lower-dimensional cell. Locally, tangential fluxes are approximated using ℝ𝕋ℕ0(K), whereas mortar fluxes and pressures using ℙ0(K).
We will not at this point place any conditions on the grid partitions, although several aspects of this will be advantageous from the perspective of computation.
6.2 Finite element spaces and the approximated problem
Let us introduce the finite element spaces necessary to write the approximated problem. We start by defining a local space for the approximated pressures, mortar fluxes, and tangential fluxes. They are given, respectively, by
where ℙ0 and ℝ𝕋ℕ0 denote the spaces of constants and lowest-order Raviart–Thomas(–Nédélec) spaces of vector functions [41, 54]. See also Fig. 4 for the degrees of freedom involved in the generic coupling between a (higher-dimensional) triangle, a mortar line segment, and a (lower-dimensional) line segment.
The composite space for the approximated mD pressure Qh ⊂ L2(Ω) and the approximated mD flux Xh ⊂ H(div; Ω, Γ) are defined respectively by
While not strictly necessary from a theoretical perspective, in the discrete setting, it is often useful to choose a finite-dimensional reconstruction operator based on the discrete spaces, and we allow for this through the notation 𝓡h,j : Λh,j → H(div; Ωi), which in practice is often further restricted to 𝓡h,j : Λh,j → Vh,ĵ. Such discrete reconstruction operators are natural for matching grids, and can also be constructed in the more general case of non-matching grids (see, e.g., [9, 16, 19]). Here Πh : Λh,j → Λ̃h,j is the L2 projection from the mortar grid on Γj to the boundary simplicial partition of Ωĵ.
We have now all the elements necessary to write the finite-dimensional approximation to the dual mixed problem (3.21).
Definition 6.1
(approximated mD dual mixed formulation). Find (𝔲h, 𝔭h) ∈ Xh × Qh such that
Due to the presence of the discrete reconstruction operator, this approximation is conforming whenever Λh,j = Λ̃h,j, i.e., for matching grids. For non-matching grids, the approximation is still convergent, subject to normal conditions on the mortar grids [19].
Remark 6.1
(conservation properties). Whenever equation (6.2b) is satisfied exactly, then equation (5.14) holds, and we have local mass conservation for matching grids. Thus, the fluxes lie in the smaller space Xh ∩ H(div; Ω, Γ;
6.3 Pressure reconstruction
Recall that Theorem 5.1 requires any approximation to the mD flux to be in H(div; Ω, Γ), whereas approximations to the mD pressure must lie in
Definition 6.2
(reconstructed pressure). We will call reconstructed pressure 𝔭̃h to any function constructed from the mD pair (𝔭h, 𝔲h) ∈ L2(Ω) × H(div; Ω, Γ) such that
Remark 6.3
(on potential reconstruction). Several techniques for obtaining 𝔭̃h are available in the literature. Arguably, the simplest option is to perform an average of the ℙ0(K) pressures on local patches and from there construct local affine ℙ1(K) functions [23]. Other techniques aim at solving first a local Neumann problem to obtain a ℙ2(K) post-processed pressure, and then apply interpolation techniques to get energy-conforming potentials [4, 5, 6, 27, 63]. Any of these choices are compatible with the bounds derived herein.
Remark 6.4
(computable estimates). Computable versions of the majorants are now readily available by setting (𝔮, 𝔳) = (𝔭̃h, 𝔲h) in (5.9), (5.13), (5.16), and (5.17).
Remark 6.5
(other locally mass-conservative methods). In addition to the MFEM scheme of the lowest-order (RT0-P0), other flux-based numerical methods such as the Mixed Virtual Element Method (MVEM) [24, 31] and Cell Centered Finite Volume Methods (CCFVM), including the Two-Point Flux Approximation (TPFA) and the Multi-Point Flux Approximation (MPFA) [1, 45], can be analyzed with our framework provided that the fluxes are interpolated in Xh} and the pressures reconstructed as indicated above. For methods without an explicit flux representation, an additional flux reconstruction step may be needed.
Remark 6.6
(superconvergence of the residual estimators). Due to Remark 6.1, the residual estimators ηR,K, LC are superconvergent for lowest-order locally mass-conservative approximations. This property is guaranteed since: (1) local Poincaré constants decay as 𝓞(hK) for simplicial elements and (2) the norm of the residual ∥fi − ∇i ⋅ vi + ∑j∈𝓢̂i νj∥K also decays as 𝓞(hK) [15]; leading to an overall rate of 𝓞(
7 Numerical validations and applications
In this section, we apply our estimators to numerical validations and benchmarks, both in two and three dimensions. To this aim, we use four different numerical methods, namely those mentioned in Remark 6.5: RT0-P0, MVEM-P0, MPFA, and TPFA. In all cases, we only consider strictly matching grids and use a low-order pressure reconstruction (recall Remark 6.3 for further discussion).
The numerical examples are implemented in the Python-based open-source software PorePy [34], using the extension package mdestimates [60], which includes the scripts of all numerical examples considered in this section.
7.1 Numerical validations
We validate the a posteriori bounds and assess their efficiency on a 1D/2D problem (Section 7.1.2) and a 2D/3D problem (Section 7.1.3), both with manufactured solutions. The geometric configuration for both problems is shown in Fig. 5. Let us denote the fracture as Ω1, the matrix as Ω2, the left interface as Γ1, and the right interface as Γ2. Further, assume the existence of an exact, smooth pressure p2(x) in Ω2. Refer to Tables 7 and 8 from Appendix D for the analytical expressions of all variables of interest.
Geometric setups used for the numerical validations. Left: A 1D fracture embedded in a 2D matrix and the exact pressure solution. Right: A 2D fracture embedded in a 3D matrix.
7.1.1 Efficiency indices
Efficiency indices are used to assess the performance of the approximations when exact solutions are available. They are defined as the ratio between the estimated and the exact errors. Here, we consider the following efficiency indices.
Definition 7.1
(efficiency indices). Let α = NC, SC, LC, EC and let 𝔭 ∈
Remark 7.1
Optimal efficiency indices (equal to 1) are obtained when the approximations match the exact solutions. Moreover, in general, the efficiency indices satisfy the bounds:
For the final term, we note that since ηR;α ⩽ 𝓜α, then for α = NC, SC the total efficiency index satisfies I𝔭,𝔲;α ⩽ 3, while for local conservation I𝔭,𝔲;LC ⩽ 2 + 𝓞(h2) and finally for exact conservation I𝔭,𝔲;EC ⩽ 2.
7.1.2 Two-dimensional validation
For our first validation, we consider the 1D/2D case as shown in the left panel of Fig. 5. This validation has two purposes: (1) compare the majorants and efficiency indices obtained using global (no mass-conservation) and local (local mass-conservation) Poincaré–Friedrichs constants, and (2) show the different errors associated with subdomains and interfaces.
To this aim, we consider four levels of successively refined combinations of mesh sizes, characterized by hcoup = h∂1Ω2 = hΓ1 = hΩ1 = hΓ2 = h∂2Ω2. The global Poincaré constant is obtained numerically by solving the associated eigenvalue problem (see, e.g., [48]), giving a value of CΩ,Γ ≈ 0.2251.
Majorants for the primal, dual, and primal–dual variables are shown in Table 1. We can see that all majorants reflect the convergence tendency of the numerical methods, and in particular (as is well-known), we identify that the TPFA approximation performs relatively poorly on this problem. As expected, the majorants obtained exploiting the local conservation properties of the methods are sharper than the ones obtained using global weights, both in absolute value and in terms of efficiency index.
Two-dimensional validation: Majorants and efficiency indices. The results for
| hcoup | I𝔭;NC | I𝔭;LC | I𝔲;NC | I𝔲;LC | I𝔭, 𝔲;NC | I𝔭,𝔲;LC | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | 0.05 | 5.86e-02 | 4.36e-02 | 1.33e-01 | 8.83e-02 | 1.46 | 1.08 | 4.09 | 3.04 | 1.89 | 1.59 |
| 0.025 | 3.01e-02 | 2.17e-02 | 6.89e-02 | 4.38e-02 | 1.49 | 1.07 | 4.18 | 3.02 | 1.91 | 1.58 | |
| 0.0125 | 1.52e-02 | 1.08e-02 | 3.48e-02 | 2.17e-02 | 1.50 | 1.07 | 4.22 | 3.00 | 1.92 | 1.57 | |
| 0.00625 | 7.65e-03 | 5.37e-03 | 1.76e-02 | 1.08e-02 | 1.52 | 1.07 | 4.25 | 2.98 | 1.93 | 1.57 | |
| MVEM-P0 | 0.05 | 6.18e-02 | 4.68e-02 | 1.40e-01 | 9.47e-02 | 1.42 | 1.07 | 4.31 | 3.26 | 1.89 | 1.60 |
| 0.025 | 3.10e-02 | 2.27e-02 | 7.08e-02 | 4.56e-02 | 1.46 | 1.07 | 4.31 | 3.15 | 1.91 | 1.59 | |
| 0.0125 | 1.54e-02 | 1.10e-02 | 3.53e-02 | 2.22e-02 | 1.49 | 1.07 | 4.29 | 3.07 | 1.92 | 1.58 | |
| 0.00625 | 7.72e-03 | 5.44e-03 | 1.77e-02 | 1.09e-02 | 1.51 | 1.06 | 4.28 | 3.02 | 1.92 | 1.57 | |
| MPFA | 0.05 | 5.91e-02 | 4.41e-02 | 1.34e-01 | 8.93e-02 | 1.46 | 1.09 | 4.12 | 3.07 | 1.89 | 1.59 |
| 0.025 | 3.03e-02 | 2.19e-02 | 6.92e-02 | 4.41e-02 | 1.49 | 1.08 | 4.20 | 3.04 | 1.91 | 1.58 | |
| 0.0125 | 1.52e-02 | 1.08e-02 | 3.49e-02 | 2.18e-02 | 1.50 | 1.07 | 4.23 | 3.01 | 1.92 | 1.57 | |
| 0.00625 | 7.66e-03 | 5.38e-03 | 1.76e-02 | 1.08e-02 | 1.52 | 1.07 | 4.25 | 2.99 | 1.93 | 1.57 | |
| TPFA | 0.05 | 6.67e-02 | 5.17e-02 | 1.50e-01 | 1.04e-01 | 1.54 | 1.19 | 3.09 | 2.39 | 1.84 | 1.58 |
| 0.025 | 3.74e-02 | 2.90e-02 | 8.35e-02 | 5.83e-02 | 1.68 | 1.31 | 2.36 | 1.83 | 1.78 | 1.52 | |
| 0.0125 | 2.64e-02 | 2.20e-02 | 5.73e-02 | 4.41e-02 | 1.82 | 1.52 | 1.64 | 1.36 | 1.63 | 1.44 | |
| 0.00625 | 1.37e-02 | 1.15e-02 | 2.98e-02 | 2.30e-02 | 1.64 | 1.37 | 1.82 | 1.52 | 1.64 | 1.44 |
Further inspection shows that efficiency indices lie within the expected bounds discussed in Remark 7.1. In particular, efficiency indices for the primal variable using local weights are very accurate, and only a ∼7% deviation with respect to the actual error (for the finest grid) is observed in the case of RT0-P0, MVEM-P0, and MPFA. For TPFA, the efficiency index is worse, as a consequence of the flux approximation being worse. Efficiency indices for the dual variable are in general larger than the ones obtained for the primal variable; this is to be expected for mixed-dual approximations with the relatively simple pressure reconstruction, where the approximated fluxes have relatively good accuracy as compared to the reconstructed pressures. Finally, efficiency indices for the primal–dual variable are less than 2 for all methods in consideration.
Considering now the local error indicators, shown in Table 2, we note that diffusive errors decrease linearly for the matrix, fracture, and interfaces. Likewise, residual errors for the matrix and fracture decrease linearly when the global Poincaré–Friedrichs constant is used. When the local Poincaré-Freidrich constants are used, the residual estimators for the matrix and the fracture decrease quadratically, which goes in agreement with the super-convergent properties discussed in Remark 6.6.
Two-dimensional validation: Subdomain and interface errors.
| hcoup | εDF, Ω2 | εR, Ω2;NC | εR, Ω2;LC | εDF, Ω1 | εR, Ω1;NC | εR, Ω1;LC | εDF, Γ1 | εDF, Γ2 | |
|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | 0.05 | 4.24e-02 | 1.41e-02 | 1.00e-03 | 2.26e-03 | 7.99e-03 | 5.65e-04 | 1.89e-04 | 1.89e-04 |
| 0.025 | 2.14e-02 | 7.73e-03 | 3.02e-04 | 1.14e-03 | 4.01e-03 | 1.42e-04 | 9.03e-05 | 9.15e-05 | |
| 0.0125 | 1.07e-02 | 4.00e-03 | 7.28e-05 | 5.70e-04 | 2.01e-03 | 3.55e-05 | 4.41e-05 | 4.41e-05 | |
| 0.00625 | 5.34e-03 | 2.07e-03 | 1.91e-05 | 2.85e-04 | 1.00e-03 | 8.87e-06 | 2.20e-05 | 2.20e-05 | |
| MVEM-P0 | 0.05 | 4.55e-02 | 1.41e-02 | 1.00e-03 | 3.25e-03 | 7.99e-03 | 5.65e-04 | 2.52e-04 | 2.52e-04 |
| 0.025 | 2.23e-02 | 7.73e-03 | 3.02e-04 | 1.32e-03 | 4.01e-03 | 1.42e-04 | 1.00e-04 | 1.03e-04 | |
| 0.0125 | 1.09e-02 | 4.00e-03 | 7.28e-05 | 5.98e-04 | 2.01e-03 | 3.55e-05 | 4.50e-05 | 4.50e-05 | |
| 0.00625 | 5.41e-03 | 2.07e-03 | 1.91e-05 | 2.89e-04 | 1.00e-03 | 8.87e-06 | 2.21e-05 | 2.21e-05 | |
| MPFA | 0.05 | 4.29e-02 | 1.41e-02 | 1.00e-03 | 2.52e-03 | 7.99e-03 | 5.65e-04 | 2.05e-04 | 2.05e-04 |
| 0.025 | 2.15e-02 | 7.73e-03 | 3.02e-04 | 1.18e-03 | 4.01e-03 | 1.42e-04 | 9.24e-05 | 9.35e-05 | |
| 0.0125 | 1.07e-02 | 4.00e-03 | 7.28e-05 | 5.77e-04 | 2.01e-03 | 3.55e-05 | 4.44e-05 | 4.44e-05 | |
| 0.00625 | 5.36e-03 | 2.07e-03 | 1.91e-05 | 2.86e-04 | 1.00e-03 | 8.87e-06 | 2.20e-05 | 2.20e-05 | |
| TPFA | 0.05 | 5.04e-02 | 1.41e-02 | 1.00e-03 | 2.52e-03 | 7.99e-03 | 5.65e-04 | 1.87e-04 | 1.89e-04 |
| 0.025 | 2.86e-02 | 7.73e-03 | 3.02e-04 | 1.18e-03 | 4.01e-03 | 1.42e-04 | 9.42e-05 | 9.23e-05 | |
| 0.0125 | 2.19e-02 | 4.00e-03 | 7.28e-05 | 5.77e-04 | 2.01e-03 | 3.55e-05 | 4.47e-05 | 4.46e-05 | |
| 0.00625 | 1.14e-02 | 2.07e-03 | 1.91e-05 | 2.86e-04 | 1.00e-03 | 8.87e-06 | 2.20e-05 | 2.21e-05 |
7.1.3 Three-dimensional validation
For our next numerical validation, we employ the 2D/3D configuration from the right panel of Fig. 5. We repeat the same analysis from the previous section, and investigate four refinement levels. The mixed-dimensional Poincaré constant for this configuration corresponds to a value of CΩ,Γ ≈ 0.1838. The results are shown in Tables 3 and 4. As in the previous validation, we can see that the majorants capture the local and global convergence tendency of all numerical methods. Again, RT0-P0, MVEM-P0, and MPFA give quite similar results, whereas TPFA showcase larger errors. As expected, efficiency indices again lie within the stated bounds from Remark 7.1.
Three-dimensional validation: Majorants and efficiency indices. The results for
| hcoup | I𝔭;NC | I𝔭;LC | I𝔲;NC | I𝔲;LC | I𝔭, 𝔲;NC | I𝔭,𝔲;LC | |||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | 0.2625 | 2.85e-01 | 2.36e-01 | 6.21e-01 | 4.73e-01 | 1.25 | 1.03 | 4.14 | 3.43 | 1.72 | 1.43 |
| 0.1720 | 1.94e-01 | 1.62e-01 | 4.20e-01 | 3.24e-01 | 1.23 | 1.03 | 4.22 | 3.53 | 1.73 | 1.50 | |
| 0.0827 | 1.07e-01 | 8.69e-02 | 2.33e-01 | 1.74e-01 | 1.28 | 1.04 | 3.61 | 2.94 | 1.70 | 1.48 | |
| 0.0418 | 5.62e-02 | 4.58e-02 | 1.23e-01 | 9.16e-02 | 1.25 | 1.02 | 3.16 | 2.58 | 1.63 | 1.43 | |
| MVEM-P0 | 0.2625 | 2.89e-01 | 2.40e-01 | 6.28e-01 | 4.80e-01 | 1.24 | 1.03 | 4.19 | 3.48 | 1.72 | 1.44 |
| 0.1720 | 1.96e-01 | 1.64e-01 | 4.24e-01 | 3.28e-01 | 1.23 | 1.03 | 4.26 | 3.57 | 1.73 | 1.50 | |
| 0.0827 | 1.08e-01 | 8.80e-02 | 2.35e-01 | 1.76e-01 | 1.27 | 1.04 | 3.65 | 2.98 | 1.70 | 1.48 | |
| 0.0418 | 5.66e-02 | 4.62e-02 | 1.24e-01 | 9.23e-02 | 1.25 | 1.02 | 3.18 | 2.60 | 1.63 | 1.44 | |
| MPFA | 0.2625 | 2.90e-01 | 2.40e-01 | 6.29e-01 | 4.82e-01 | 1.25 | 1.03 | 4.08 | 3.38 | 1.72 | 1.43 |
| 0.1720 | 1.98e-01 | 1.66e-01 | 4.28e-01 | 3.32e-01 | 1.23 | 1.03 | 4.22 | 3.54 | 1.73 | 1.50 | |
| 0.0827 | 1.09e-01 | 8.90e-02 | 2.37e-01 | 1.78e-01 | 1.27 | 1.04 | 3.64 | 2.98 | 1.70 | 1.49 | |
| 0.0418 | 5.69e-02 | 4.65e-02 | 1.24e-01 | 9.30e-02 | 1.25 | 1.02 | 3.18 | 2.60 | 1.63 | 1.44 | |
| TPFA | 0.2625 | 3.84e-01 | 3.35e-01 | 8.17e-01 | 6.70e-01 | 1.24 | 1.08 | 2.13 | 1.86 | 1.48 | 1.28 |
| 0.1720 | 2.95e-01 | 2.63e-01 | 6.23e-01 | 5.27e-01 | 1.38 | 1.23 | 1.66 | 1.48 | 1.44 | 1.30 | |
| 0.0827 | 2.22e-01 | 2.02e-01 | 4.63e-01 | 4.04e-01 | 1.62 | 1.48 | 1.40 | 1.28 | 1.45 | 1.35 | |
| 0.0418 | 2.08e-01 | 1.97e-01 | 4.26e-01 | 3.95e-01 | 1.76 | 1.67 | 1.29 | 1.23 | 1.46 | 1.41 |
Three-dimensional validation: Subdomain and interface errors.
| hcoup | εDF, Ω2 | εR, Ω2;NC | εR, Ω2;LC | εDF, Ω1 | εR, Ω1;NC | εR, Ω1;LC | εDF, Γ1 | εDF, Γ2 | |
|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | 0.2625 | 2.35e-01 | 4.73e-02 | 2.55e-02 | 4.67e-03 | 1.56e-02 | 6.70e-03 | 5.05e-03 | 5.00e-03 |
| 0.1720 | 1.62e-01 | 3.08e-02 | 1.10e-02 | 4.53e-03 | 9.01e-03 | 2.04e-03 | 1.37e-03 | 1.37e-03 | |
| 0.0827 | 8.69e-02 | 1.91e-02 | 3.91e-03 | 2.72e-03 | 5.17e-03 | 6.70e-04 | 4.07e-04 | 4.09e-04 | |
| 0.0418 | 4.58e-02 | 1.01e-02 | 1.06e-03 | 1.40e-03 | 2.64e-03 | 1.70e-04 | 1.08e-04 | 1.09e-04 | |
| MVEM-P0 | 0.2625 | 2.39e-01 | 4.73e-02 | 2.55e-02 | 5.78e-03 | 1.56e-02 | 6.70e-03 | 5.54e-03 | 5.49e-03 |
| 0.1720 | 1.64e-01 | 3.08e-02 | 1.10e-02 | 5.56e-03 | 9.01e-03 | 2.04e-03 | 1.50e-03 | 1.50e-03 | |
| 0.0827 | 8.79e-02 | 1.91e-02 | 3.91e-03 | 3.05e-03 | 5.17e-03 | 6.70e-04 | 4.40e-04 | 4.41e-04 | |
| 0.0418 | 4.61e-02 | 1.01e-02 | 1.06e-03 | 1.46e-03 | 2.64e-03 | 1.70e-04 | 1.13e-04 | 1.14e-04 | |
| MPFA | 0.2625 | 2.40e-01 | 4.73e-02 | 2.55e-02 | 5.04e-03 | 1.56e-02 | 6.70e-03 | 5.93e-03 | 5.86e-03 |
| 0.1720 | 1.66e-01 | 3.08e-02 | 1.10e-02 | 4.86e-03 | 9.01e-03 | 2.04e-03 | 1.56e-03 | 1.55e-03 | |
| 0.0827 | 8.89e-02 | 1.91e-02 | 3.91e-03 | 2.82e-03 | 5.17e-03 | 6.70e-04 | 4.61e-04 | 4.62e-04 | |
| 0.0418 | 4.65e-02 | 1.01e-02 | 1.06e-03 | 1.41e-03 | 2.64e-03 | 1.70e-04 | 1.14e-04 | 1.16e-04 | |
| TPFA | 0.2625 | 3.34e-01 | 4.73e-02 | 2.55e-02 | 4.88e-03 | 1.56e-02 | 6.70e-03 | 6.04e-03 | 5.11e-03 |
| 0.1720 | 2.63e-01 | 3.08e-02 | 1.10e-02 | 4.86e-03 | 9.01e-03 | 2.04e-03 | 1.35e-03 | 1.29e-03 | |
| 0.0827 | 2.02e-01 | 1.91e-02 | 3.91e-03 | 2.85e-03 | 5.17e-03 | 6.70e-04 | 4.50e-04 | 4.39e-04 | |
| 0.0418 | 1.97e-01 | 1.01e-02 | 1.06e-03 | 1.46e-03 | 2.64e-03 | 1.70e-04 | 1.02e-04 | 1.02e-04 |
7.2 Numerical applications
We now apply our estimators to numerical approximations of challenging problems solving the equations of incompressible flow in fractured porous media. Importantly, since source terms are zero in both applications, by applying matching grids the residual errors are zero, and we are in the setting of having an exact conservation property from the numerical approximation. From Remark 7.1, we then know that the efficiency index for the primal–dual error will be less than 2; even if the exact solution and error are both unknown.
7.2.1 Two-dimensional application
In this numerical experiment, we consider the benchmark case 3b from [29]. As shown in the left panel of Fig. 1, the domain consists of ten (partially intersecting) fractures embedded in a unit square matrix. The exact fracture coordinates can be found in Appendix C of [29]. Fractures 4 and 5 represent blocking fractures (𝓚 = 10−4 and κ = 1) whereas the others represent conductive fractures (𝓚 = 104 and κ = 108). The matrix permeability is set to one. A linear pressure drop is imposed from left (p = 4) to right (p = 1), whereas no flux is prescribed at the top and bottom of the domain.
The benchmark establishes three refinement levels; coarse, intermediate, and fine, with approximately 1500, 4200, and 16000 two-dimensional cells. The structure of the local contributions to the majorant (confer, e.g., equation (5.7)) are shown in Fig. 6, based on the approximate solution obtained by the MPFA discretization.
In Table 5, we show the errors bounds for the three refinement levels. To avoid numbering domains and interfaces, we refer to the matrix error as εΩ2, EC, and group the fracture and interface errors by conductive and blocking. For example, εΩ1, C, EC refers to the sum of the errors of 1D conductive fractures.
Error estimates for the two-dimensional application. The results for
| Mesh | εΩ2;EC | εΩ1;EC, C | εΩ1;EC, B | εΓ1, C | εΓ1, B | εΓ0 | |||
|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | Coarse | 7.39e-01 | 2.93e-01 | 2.98e-04 | 3.13e+03 | 1.52e-01 | 2.24e+01 | 9.94e+02 | 1.99e+03 |
| Intermediate | 5.95e-01 | 1.90e-01 | 2.77e-04 | 1.95e+03 | 1.00e-01 | 1.79e+01 | 6.20e+02 | 1.24e+03 | |
| Fine | 4.30e-01 | 1.07e-01 | 2.78e-04 | 9.79e+02 | 5.15e-02 | 1.22e+01 | 3.15e+02 | 6.30e+02 | |
| MVEM-P0 | Coarse | 7.29e-01 | 3.51e-01 | 1.44e-04 | 3.10e+03 | 1.46e-01 | 4.41e+01 | 9.84e+02 | 1.97e+03 |
| Intermediate | 5.91e-01 | 2.23e-01 | 1.27e-04 | 1.94e+03 | 9.43e-02 | 3.14e+01 | 6.17e+02 | 1.23e+03 | |
| Fine | 4.28e-01 | 1.24e-01 | 1.18e-04 | 9.78e+02 | 4.80e-02 | 2.02e+01 | 3.15e+02 | 6.29e+02 | |
| MPFA | Coarse | 7.39e-01 | 3.13e-01 | 1.72e-04 | 3.03e+03 | 1.43e-01 | 3.39e+01 | 9.63e+02 | 1.93e+03 |
| Intermediate | 5.98e-01 | 2.01e-01 | 1.54e-04 | 1.89e+03 | 9.18e-02 | 2.55e+01 | 6.00e+02 | 1.20e+03 | |
| Fine | 4.33e-01 | 1.12e-01 | 1.46e-04 | 9.49e+02 | 4.71e-02 | 1.68e+01 | 3.05e+02 | 6.10e+02 | |
| TPFA | Coarse | 7.52e-01 | 3.05e-01 | 1.76e-04 | 3.19e+03 | 1.48e-01 | 3.67e+01 | 1.01e03 | 2.02e+03 |
| Intermediate | 6.08e-01 | 1.96e-01 | 1.51e-04 | 1.95e+03 | 9.41e-02 | 2.61e+01 | 6.12e+02 | 1.22e+03 | |
| Fine | 4.45e-01 | 1.09e-01 | 1.60e-04 | 1.00e+03 | 4.84e-02 | 1.86e+01 | 3.23e+02 | 6.46e+02 |
An important observation is that the persistent reduction of the majorant
The error estimates suggest that the contribution to the overall error bounds are concentrated, primarily, on highly conductive interfaces (see the column corresponding to εΓ1,C). On a more qualitative note, Figure 6 suggests that subdomain diffusive errors are concentrated at the fracture tips and fracture intersections, which is where singularities may typically be encountered [19].
Two-dimensional benchmark problem and the errors associated with the matrix and fractures for the coarse (left), intermediate (center), and fine (right) grid resolutions. Fractures 4 and 5 are blocking, whereas the others are conductive. The local bounds were obtained using MPFA. The results suggest that subdomain diffusive errors are concentrated around fracture tips and fracture intersections.
7.2.2 Three-dimensional application
Our last numerical application is based on a modified version of the three-dimensional benchmark 2.1 from [14]. The domain consists of nine intersecting fractures embedded in a unit cube, as shown in the middle panel of Fig. 1. This results in an intricate network with 106 subdomains and 270 interfaces of different dimensionality.
The original benchmark imposes an inlet flux (purple lower corner {u} = -1) and an outlet pressure (pink upper corner p = 1), and for the rest of the external boundaries null flux. Since we have only detailed our results for zero Neumann boundary conditions, we have replaced the inlet flux by a constant pressure condition (p = 1) and modified the value of the outlet pressure (p = 0). The benchmark assigns heterogeneous permeability to the matrix subdomain, whereas the fractures are assumed to be highly conductive. For the complete description of the benchmark, we refer to [14], and for an impression on how the contributions to the majorant are distributed (see Fig. 7). Here we show the error estimates for the whole fracture network obtained with RT0-P0, where it becomes evident that the subdomain diffusive errors are concentrated at the inlet and outlet boundaries; refinement efforts should therefore focus on these regions.
As in Section 7.2.1, we collect the local errors of subdomains and interfaces of equal dimensionality. The results are summarized in Table 6. As in the previous cases, we have local and global convergence for all four numerical methods. Again, RT0-P0, MVEM-P0, and MPFA show very similar results, while TPFA show larger errors.
Error estimates for the three-dimensional application. The results for
| Mesh | εΩ3;EC | εΩ2;EC | εΩ1;EC | εΓ2 | εΓ1 | εΓ0 | |||
|---|---|---|---|---|---|---|---|---|---|
| RT0-P0 | Coarse | 6.17e-01 | 5.81e-04 | 3.16e-04 | 9.87e+02 | 3.63e-02 | 3.31e-02 | 5.03e+02 | 1.01e+03 |
| Intermediate | 4.55e-01 | 4.61e-04 | 1.58e-04 | 7.75e+01 | 8.86e-03 | 8.35e-04 | 3.40e+01 | 6.81e+01 | |
| Fine | 3.86e-01 | 2.55e-04 | 9.60e-05 | 2.26e+01 | 4.63e-03 | 4.34e-04 | 1.07e+01 | 2.14e+01 | |
| MVEM-P0 | Coarse | 6.07e-01 | 6.99e-04 | 2.77e-04 | 9.54e+02 | 7.48e-02 | 6.38e-02 | 4.66e+02 | 9.33e+02 |
| Intermediate | 4.55e-01 | 4.63e-04 | 1.65e-04 | 8.19e+01 | 9.96e-03 | 4.59e-03 | 3.59e+01 | 7.18e+01 | |
| Fine | 3.86e-01 | 2.46e-04 | 9.17e-05 | 2.33e+01 | 4.00e-03 | 1.75e-03 | 1.11e+01 | 2.22e+01 | |
| MPFA | Coarse | 6.07e-01 | 7.00e-04 | 3.15e-04 | 1.05e+03 | 4.61e-02 | 1.69e-02 | 5.24e+02 | 1.05e+03 |
| Intermediate | 4.46e-01 | 4.88e-04 | 1.61e-04 | 8.42e+01 | 7.72e-03 | 2.31e-03 | 3.71e+01 | 7.42e+01 | |
| Fine | 3.77e-01 | 2.53e-04 | 9.04e-05 | 2.37e+01 | 2.82e-03 | 9.36e-04 | 1.12e+01 | 2.24e+01 | |
| TPFA | Coarse | 6.32e-01 | 4.72e-04 | 2.26e-04 | 7.92e+02 | 4.21e-02 | 1.34e-02 | 3.76e+02 | 7.52e+02 |
| Intermediate | 4.48e-01 | 6.27e-04 | 1.40e-04 | 1.47e+02 | 1.56e-02 | 2.32e-03 | 6.82e+01 | 1.36e+02 | |
| Fine | 4.07e-01 | 5.82e-04 | 8.72e-05 | 4.60e+01 | 7.97e-03 | 1.05e-03 | 2.04e+01 | 4.08e+01 |
As in the 2D case discussed above, the persistent reduction of the majorant
Subdomain diffusive error contributions to the majorant for the fine grid resolution obtained with RT0-P0.
8 Conclusion
In this paper, we obtained a posteriori error estimates for mixed-dimensional elliptic equations. Depending upon the level of accuracy at which residual balances can be approximated, we have derived four concrete versions of the majorant; i.e., for no mass-conservative, subdomain mass-conservative, locally mass-conservative, and point-wise mass-conservative approximations. Furthermore, we have demonstrated both theoretically and numerically that sharper bounds can be obtained (for locally mass-conservative methods) using local Poincaré constants instead of the global ones.
Our bounds have been thoroughly tested with numerical approximations obtained with four locally mass-conservative methods of the lowest-order, namely: RT0-P0, MVEM-P0, MPFA, and TPFA. We performed a detailed efficiency analysis comparing the use of global and local Poincaré–Friedrichs constants in two and three dimensions. In both validations, our upper bounds reflected the optimal convergence rates of the numerical methods. In addition, we applied our bounds to two- and three-dimensional community benchmark problems exhibiting challenging fracture networks. Again, in both cases, the bounds reflected the limitations and the convergence rates of the methods satisfactory.
To the best of our knowledge, the bounds obtained here are the first of their kind to provide a practical tool to measure the error in numerical approximations to the equations modeling the incompressible, single-phase flow in generic fractured porous media.
-
Funding: J. Varela was funded by VISTA — a basic research program in collaboration between The Norwegian Academy of Science and Letters, and Equinor. Additionally, this work was supported in part through the Norwegian Research Council grant 250223. The authors would like to thank W. M. Boon for helpful discussions on this topic.
References
[1] I. Aavatsmark, An introduction to multipoint flux approximations for quadrilateral grids, Comput. Geosci. 6 (2002), No. 3-4, 405–432.Search in Google Scholar
[2] R. A. Adams and J. J. F. Fournier, Sobolev Spaces, 2nd ed., Pure and Applied Mathematics, Vol. 140, Elsevier/Academic Press, Amsterdam, 2003.Search in Google Scholar
[3] E. Ahmed, J. Jaffr’, and J. E. Roberts, A reduced fracture model for two-phase flow with different rock types, Math. Comput. Simulation 137 (2017), 49–70.10.1016/j.matcom.2016.10.005Search in Google Scholar
[4] E. Ahmed, J. M. Nordbotten, and F. A. Radu, Adaptive asynchronous time-stepping, stopping criteria, and a posteriori error estimates for fixed-stress iterative schemes for coupled poromechanics problems, J. Comput. Appl. Math. 364 (2020), 112312, 25.10.1016/j.cma.2018.12.016Search in Google Scholar
[5] E. Ahmed, F. A. Radu, and J. M. Nordbotten, Adaptive poromechanics computations based on a posteriori error estimates for fully mixed formulations of Biot’s consolidation model, Comput. Methods Appl. Mech. Engrg. 347 (2019), 264–294.10.1016/j.cam.2019.06.028Search in Google Scholar
[6] M. Ainsworth, A posteriori error estimation for lowest order Raviart–Thomas mixed finite elements, SIAM J. Sci. Comput. 30 (2007/08), No. 1, 189–204.10.1137/06067331XSearch in Google Scholar
[7] C. Alboin, J. Jaffr’, J. E. Roberts, and C. Serres, Modeling fractures as interfaces for 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.10.1090/conm/295/04999Search in Google Scholar
[8] S. S. Antman, Nonlinear Problems of Elasticity, Applied Mathematical Sciences, Vol. 107, Springer-Verlag, New York, 1995.10.1007/978-1-4757-4147-6_1Search in Google Scholar
[9] T. Arbogast, L. C. Cowsar, M. F. Wheeler, and I. Yotov, Mixed finite element methods on nonmatching multiblock grids, SIAM J. Numer. Analysis 37 (2000), No. 4, 1295–1315.10.1137/S0036142996308447Search in Google Scholar
[10] D. N. Arnold, Finite Element Exterior Calculus, CBMS-NSF Regional Conference Series in Applied Mathematics, Vol. 93, SIAM, Philadelphia, PA, 2018.10.1137/1.9781611975543.ch1Search in Google Scholar
[11] I. Babuška and W. C. Rheinboldt, Error estimates for adaptive finite element computations, SIAM J. Numer. Anal. 15 (1978), No. 4, 736–754.10.1137/0715049Search in Google Scholar
[12] M. Bebendorf, A note on the Poincar’ inequality for convex domains, Z. Anal. Anwendungen 22 (2003), No. 4, 751–756.10.4171/ZAA/1170Search in Google Scholar
[13] Z. Belhachmi, A posteriori error estimates for the 3D stabilized mortar finite element method applied to the Laplace equation, M2AN Math. Model. Numer. Anal. 37 (2003), No. 6, 991–1011.10.1051/m2an:2003064Search in Google Scholar
[14] I. Berre, W. M. Boon, B. Flemisch, A. Fumagalli, D. Gläser, E. Keilegavlen, A. Scotti, I. Stefansson, A. Tatomir, K. Brenner, S. Burbulla, P. Devloo, O. Duran, M. Favino, J. Hennicker, I-H. Lee, K. Lipnikov, R. Masson, K. Mosthaf, M. G. C. Nestola, C-F. Ni, K. Nikitin, P. Schädle, D. Svyatskiy, R. Yanbarisov, and P. Zulian, Verification benchmarks for single-phase flow in three-dimensional fractured porous media, Advances in Water Resources 147 (2021), 103759.10.1016/j.advwatres.2020.103759Search in Google Scholar
[15] D. Boffi, F. Brezzi, and M. Fortin, Mixed Finite Element Methods and Applications, Springer Series in Computational Mathematics, Vol. 44, Springer, Heidelberg, 2013.10.1007/978-3-642-36519-5Search in Google Scholar
[16] W. M. Boon, D. Gläser, R. Helmig, and I. Yotov, Flux-mortar mixed finite element methods on non-matching grids, SIAM J. Numer. Anal. 60 (2022), No. 3, 1193–1225. 2020.10.1007/s10596-020-10013-2Search in Google Scholar
[17] W. M. Boon and J. M. Nordbotten, Stable mixed finite elements for linear elasticity with thin inclusions, Computational Geosciences (2021), 603–620.10.1137/17M1139102Search in Google Scholar
[18] W. M. Boon, J. M Nordbotten, and J. E. Vatne, Functional analysis and exterior calculus on mixed-dimensional geometries, Annali di Matematica Pura ed Applicata (1923-) 200 (2021), No. 2, 757–789.10.1007/s10231-020-01013-1Search in Google Scholar
[19] W. M. Boon, J. M. Nordbotten, and I. Yotov, Robust discretization of flow in fractured porous media, SIAM J. Numer. Anal. 56 (2018), No. 4, 2203–2233.10.1007/s10231-020-01013-1Search in Google Scholar
[20] C. Carstensen and S. A. Funken, Constants in Clément-interpolation error and residual based a posteriori estimates in finite element methods, East-West J. Numer. Math 8 (2000), No. 3, 153–175.Search in Google Scholar
[21] H. Chen and S. Sun, A residual-based a posteriori error estimator for single-phase Darcy flow in fractured porous media, Numer. Math. 136 (2017), No. 3, 805–839.10.1007/s00211-016-0851-9Search in Google Scholar
[22] P. G. Ciarlet, Mathematical Elasticity, Vol. II, Studies in Mathematics and its Applications, Vol. 27, North-Holland Publishing Co., Amsterdam, 1997.Search in Google Scholar
[23] S. Cochez-Dhondt, S. Nicaise, and S. I. Repin, A posteriori error estimates for finite volume approximations, Math. Model. Nat. Phenom. 4 (2009), No. 1, 106–122.10.1051/mmnp/20094105Search in Google Scholar
[24] L. B. da Veiga, F. Brezzi, L. D. Marini, and A. Russo, Mixed virtual element methods for general second order elliptic problems on polygonal meshes, ESAIM Math. Model. Numer. Anal. 50 (2016), No. 3, 727–747.10.1051/m2an/2015067Search in Google Scholar
[25] C. D’Angelo and A. Quarteroni, On the coupling of 1D and 3D diffusion–reaction equations. Application to tissue perfusion problems, Math. Models Methods Appl. Sci. 18 (2008), No. 8, 1481–1504.10.1142/S0218202508003108Search in Google Scholar
[26] A. Ern, I. Smears, and M. Vohralík, Guaranteed, locally space-time efficient, and polynomial-degree robust a posteriori error estimates for high-order discretizations of parabolic problems, SIAM J. Numer. Anal. 55 (2017), No. 6, 2811–2834.10.1137/090759008Search in Google Scholar
[27] A. Ern and M. Vohralík, A posteriori error estimation based on potential and flux reconstruction for the heat equation, SIAM J. Numer. Anal. 48 (2010), No. 1, 198–223.10.1137/130950100Search in Google Scholar
[28] A. Ern and M. Vohralík, Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations, SIAM J. Numer. Anal. 53 (2015), No. 2, 1058–1081.10.1137/16M1097626Search in Google Scholar
[29] B. Flemisch, I. Berre, W. Boon, A. Fumagalli, N. Schwenck, A. Scotti, I. Stefansson, and A. Tatomir, Benchmarks for single-phase flow in fractured porous media, Advances in Water Resources 111 (2018), 239–258.10.1016/j.advwatres.2017.10.036Search in Google Scholar
[30] L. Formaggia, A. Fumagalli, A. Scotti, and P. Ruffo, A reduced model for Darcy’s problem in networks of fractures, ESAIM Math. Model. Numer. Anal. 48 (2014), No. 4, 1089–1116.10.1051/m2an/2013132Search in Google Scholar
[31] A. Fumagalli and E. Keilegavlen, Dual virtual element methods for discrete fracture matrix models, Oil & Gas Science and Technology–Revue d’IFP Energies Nouvelles 74 (2019), 41.10.2516/ogst/2019008Search in Google Scholar
[32] F. Hecht, Z. Mghazli, I. Naji, and J. E. Roberts, A residual a posteriori error estimators for a model for flow in porous media with fractures, J. Sci. Comput. 79 (2019), No. 2, 935–968.10.1007/s10915-018-0875-7Search in Google Scholar
[33] E. Hodneland, E. Hanson, O. Sævareid, G. Nævdal, A. Lundervold, V. Šoltészová, A. Z. Munthe-Kaas, A. Deistung, J. R. Reichenbach, and J. M. Nordbotten, A new framework for assessing subject-specific whole brain circulation and perfusion using MRI-based measurements and a multi-scale continuous flow model, PLoS Computational Biology 15 (2019), No. 6, e1007073.10.1371/journal.pcbi.1007073Search in Google Scholar PubMed PubMed Central
[34] E. Keilegavlen, R. Berge, A. Fumagalli, M. Starnoni, I. Stefansson, J. Varela, and I. Berre, Porepy: An open-source software for simulation of multiphysics processes in fractured porous media, Computational Geosciences 25 (2021), No. 1, 243–265.10.1007/s10596-020-10002-5Search in Google Scholar
[35] D. W. Kelly, J. P. de S. R. Gago, O. C. Zienkiewicz, and I. Babuška, A posteriori error analysis and adaptive processes in the finite element method, I. Error analysis, Internat. J. Numer. Methods Engrg. 19 (1983), No. 11, 1593–1619.10.1002/nme.1620191103Search in Google Scholar
[36] T. Koch, K. Heck, N. Schröder, H. Class, and R. Helmig, A new simulation framework for soil–root interaction, evaporation, root growth, and solute transport, Vadose Zone Journal 17 (2018), No. 1, 1–21.10.2136/vzj2017.12.0210Search in Google Scholar
[37] T. Köppl, E. Vidotto, and B. Wohlmuth, A local error estimate for the Poisson equation with a line source term. In: Numerical Mathematics and Advanced Applications—ENUMATH 2015, Lect. Notes Comput. Sci. Engrg., Vol. 112, Springer, Cham, 2016, pp. 421–429.10.1007/978-3-319-39929-4_40Search in Google Scholar
[38] S. Kurz, D. Pauly, D. Praetorius, S. I. Repin, and D. Sebastian, Functional a posteriori error estimates for boundary element methods, Numer. Math. 147 (2021), No. 4, 937–966.10.1007/s00211-021-01188-6Search in Google Scholar
[39] V. Martin, J. Jaffr’, and J. E. Roberts, Modeling fractures and barriers as interfaces for flow in porous media, SIAM J. Sci. Comput. 26 (2005), No. 5, 1667–1691.10.1137/S1064827503429363Search in Google Scholar
[40] Z. Mghazli and I. Naji, Guaranteed a posteriori error estimates for a fractured porous medium, Math. Comput. Simulation 164 (2019), 163–179.10.1016/j.matcom.2019.02.002Search in Google Scholar
[41] J.-C. N’d’lec, Mixed finite elements in R3, Numer. Math. 35 (1980), No. 3, 315–341.10.1007/BF01396415Search in Google Scholar
[42] P. Neittaanmäki and S. I. Repin, Reliable methods for computer simulation, Studies in Mathematics and its Applications, Vol. 33, Elsevier Science B.V, Amsterdam, 2004.Search in Google Scholar
[43] J. M. Nordbotten, Mixed-dimensional models for real-world applications, Snapshots of Modern Mathematics from Oberwolfach (2019), 11.Search in Google Scholar
[44] J. M. Nordbotten, W. M. Boon, A. Fumagalli, and E. Keilegavlen, Unified approach to discretization of flow in fractured porous media, Comput. Geosci. 23 (2019), No. 2, 225–237.10.1007/978-3-030-69363-3_4Search in Google Scholar
[45] J. M. Nordbotten and E. Keilegavlen, An introduction to multi-point flux (MPFA) and stress (MPSA) finite volume methods for thermo-poroelasticity, In: Polyhedral Methods in Geosciences, Springer, 2021, pp. 119–158.10.1007/s10596-018-9778-9Search in Google Scholar
[46] J. T. Oden and S. Prudhomme, Goal-oriented error estimation and adaptivity for the finite element method, Comput. Math. Appl. 41 (2001), No. 5-6, 735–756.10.1016/S0898-1221(00)00317-5Search in Google Scholar
[47] D. Pauly, Solution theory, variational formulations, and functional a posteriori error estimates for general first order systems with applications to electro-magneto-statics and more, Numer. Funct. Anal. Optim. 41 (2020), No. 1, 16–112.10.1080/01630563.2018.1490756Search in Google Scholar
[48] D. Pauly and J. Valdman, Poincar’–Friedrichs type constants for operators involving grad, curl, and div: theory and numerical experiments, Comput. Math. Appl. 79 (2020), No. 11, 3027–3067.10.1016/j.camwa.2020.01.004Search in Google Scholar
[49] L. E. Payne and H. F. Weinberger, An optimal Poincar’ inequality for convex domains, Arch. Rational Mech. Anal. 5 (1960), 286–292.10.1007/BF00252910Search in Google Scholar
[50] C. Pechstein and R. Scheichl, Weighted Poincar’ inequalities, IMA J. Numer. Anal. 33 (2013), No. 2, 652–686.10.1093/imanum/drs017Search in Google Scholar
[51] G. K. Pedersen, Analysis Now, Graduate Texts in Mathematics, Vol. 118, Springer-Verlag, New York, 1989.10.1007/978-1-4612-1007-8Search in Google Scholar
[52] G. V. Pencheva, M. Vohralík, M. F. Wheeler, and T. Wildey, Robust a posteriori error control and adaptivity for multiscale, multinumerics, and mortar coupling, SIAM J. Numer. Anal. 51 (2013), No. 1, 526–554.10.1137/110839047Search in Google Scholar
[53] M. Rathmair, On how Poincar’ inequalities imply weighted ones, Monatsh. Math. 188 (2019), No. 4, 753–763.10.1007/s00605-019-01266-wSearch in Google Scholar
[54] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, In: Mathematical aspects of finite element methods (Proc. Conf., Consiglio Naz. delle Ricerche (C.N.R.), Rome, 1975), Lecture Notes in Math., Vol. 606, 1977, pp. 292–315.10.1007/BFb0064470Search in Google Scholar
[55] S. I. Repin, A posteriori error estimation for variational problems with uniformly convex functionals, Math. Comp. 69 (2000), No. 230, 481–500.10.1090/S0025-5718-99-01190-4Search in Google Scholar
[56] S. I. Repin, Two-sided estimates of deviation from exact solutions of uniformly elliptic equations, In: Proceedings of the St. Petersburg Mathematical Society, Vol. IX, Amer. Math. Soc. Transl. Series 2, 209, Amer. Math. Soc., Providence, RI, 2003, pp. 143–171,10.1090/trans2/209/06Search in Google Scholar
[57] S. I. Repin, A Posteriori Estimates for Partial Differential Equations, Radon Series on Computational and Applied Mathematics, Vol. 4, Walter de Gruyter GmbH & Co. KG, Berlin, 2008.10.1515/9783110203042Search in Google Scholar
[58] S. I. Repin, Computable majorants of constants in the Poincar’ and Friedrichs inequalities, J. Math. Sci. (N.Y.) 186 (2012), No. 2, 307–321.10.1007/s10958-012-0987-9Search in Google Scholar
[59] S. I. Repin, S. Sauter, and A. Smolianski, Two-sided a posteriori error estimates for mixed formulations of elliptic problems, SIAM J. Numer. Anal. 45 (2007), No. 3, 928–945.10.1137/050641533Search in Google Scholar
[60] J. Varela, jhabriel/mixdim-estimates: v1.4, March 2022.Search in Google Scholar
[61] R. Verfürth, A review of a posteriori error estimation techniques for elasticity problems, Comput. Methods Appl. Mech. Engrg. 176 (1999), No. 1-4, 419–440.10.1016/S0045-7825(98)00347-8Search in Google Scholar
[62] M. Vohralík, On the discrete Poincar’–Friedrichs inequalities for nonconforming approximations of the Sobolev space H1, Numer. Funct. Anal. Optim. 26 (2005), No. 7-8, 925–952.10.1080/01630560500444533Search in Google Scholar
[63] M. Vohralík, Unified primal formulation-based a priori and a posteriori error analysis of mixed finite element methods, Math. Comp. 79 (2010), No. 272, 2001–2032.10.1090/S0025-5718-2010-02375-0Search in Google Scholar
[64] M. F. Wheeler and I. Yotov, A posteriori error estimates for the mortar mixed finite element method, SIAM J. Numer. Anal. 43 (2005), No. 3, 1021–1042.10.1137/S0036142903431687Search in Google Scholar
[65] B. I. Wohlmuth, Hierarchical a posteriori error estimators for mortar finite element methods with Lagrange multipliers, SIAM J. Numer. Anal. 36 (1999), No. 5, 1636–1658.10.1137/S0036142997330512Search in Google Scholar
[66] O. C. Zienkiewicz and J. Z. Zhu, A simple error estimator and adaptive procedure for practical engineering analysis, Int. J. Numer. Methods Engrg. 24 (1987), No. 2, 337–357.10.1002/nme.1620240206Search in Google Scholar
A Derivation of variational formulations
Here, we present the derivations for the primal and dual variational formulations for the case of a single fracture immersed in a matrix.
A.1 Derivation of the primal weak form for a single fracture
Substitute (2.1b) into (2.1a), multiply each term by q2 ∈
Using integration by parts, the first term of (A.1) can be expressed as
Here, we use the internal boundary conditions (2.1c) and (2.1d) and the definition of the mortar fluxes (2.3a) and (2.3b). Analogously, integration by parts allows us to write the second term of (A.1) as
Note that the boundary terms vanish due to the choice of boundary conditions.
Finally, we note that the third and fourth terms from (A.1) can be equivalently written as
The proof is completed by substituting (A.2), (A.3), and (A.3) into (A.1) and grouping common terms.
A.2 Derivation of the dual weak form for a single fracture
Let us start with (2.13a). Multiply respectively (2.1b) and (2.2b) by v0,2 ∈ V0,2 and v0, 1 ∈ V0,1, integrate over the subdomains Ω2 and Ω1, use integration by parts to obtain
Adding together (A.5) and (A.6) gives (2.13a). We now focus on (2.13b). First, we use (2.1b) and multiply by the test functions 𝓡jνj with νj ∈ L2(Γj) for j ∈ {1, 2}, integrate over Ω2, and apply integration by parts, to obtain
Next, we multiply the interface laws (2.3a) and (2.3b) by ν1 and ν2, respectively, to get
After adding (A.7) and (A.8) and canceling common terms, we obtain (2.13b). Finally, to obtain (2.13c), we multiply (2.1a) by q2 ∈ L2(Ω2) and (2.2a) by q1 ∈ L2(Ω1), and integrate over their respective subdomains, and add the resulting equations.
B Proof of Theorem 2.1
Here, we present the proof of the upper bound of the error for the primal variable, for the case of a single fracture immersed in a matrix.
Proof
Start by computing the difference between p = [p1, p2] ∈
By noticing that the first three terms of (B.1) add up to the right-hand side of (2.6), and adding the identity
valid for any v0 ∈ V0 and ν ∈ L2(Γ) to (B.1), we obtain
Recognizing that since 𝓚2 is symmetric positive definite, it can be expressed as
Finally, applying the permeability-weighted Poincaré–Friedrichs inequality (4.1a) to the terms ∥p1 − q1∥Ω1 and ∥p2 − q2∥Ω2, the proof of the theorem is completed.□
C Proof of Theorem 5.1
Here, we present the proof of our main theorem, which deals with the general abstract estimates in a mixed-dimensional setting.
Proof
(1) The proof for the bounds for the mD primal variable follows the one presented in Appendix B, modulo its generalization to the mD setting and the use of weighted norms on the residual terms. Start by computing the difference between any
Here, we used (4.10), (5.3), and added the fact that 𝔇⋅ and 𝔻 are adjoints.
By exploiting the orthogonality property (4.7) and then introducing the weights to the second and third terms, (C.1) can be equivalently written as:
Finally, applying the Cauchy–Schwarz inequality to the first and second terms of (C.2), and then the norm definitions (4.10), (4.11), and (4.8), we arrive at the desired bound:
(2) The proof for the bounds for the dual variable is given next. We remark that an alternative proof based on a generalized abstract estimate (see [63], Theorem 6.1) can be used to obtain equivalent upper bounds after its generalization to the mD setting.
We start by adding the square of the primal and dual error to obtain:
Here, we used the norm definitions (4.10) and (4.11) together with the mD constitutive relationship (3.12a).
Using partial integration, mass conservation (3.12b), and the orthogonality property (4.7), the second term of (C.4) can be equivalently written as
Using the Cauchy–Schwarz inequality twice and the definition of the weighted norms (4.8), (C.5) can be estimated as
Substituting (C.6) into (C.4) and applying the Cauchy–Schwarz inequality to the first term, we arrive at
from which we conclude that (5.4) indeed holds.
(3) To prove the upper bound for the primal–dual pair, we choose an arbitrary pair
where we use the bounds (5.3) and (5.4).
For the proof of the lower bound, we start from the definition of the majorant, to get
This completes the proof for the two-sided bounds and the abstract theorem.□
D Exact solutions to numerical validations
Herein, we provide the exact expressions for the pressure, velocities, mortar fluxes, and source terms for the numerical validations presented in Section 7.1.
We will conveniently define the following quantities for notational compactness:
where x = [x1, x2, x3].
D.1 Exact solutions for the 1D/2D validation
The matrix subdomain Ω2 is decomposed into three regions, i.e.,
Let us now define the distance function d(x) from Ω2 to Ω1. That is,
and the bubble function ω(x):
In Table 7, we include the exact solutions for all the variables of interest. Note that the parameter n controls the regularity of the solution. For this particular validation, a value of n = 1.5 was adopted.
Exact solutions for the 1D/2D validation.
| p2 = | dn+1 + ω d | |
|---|---|---|
| Ω2 ∖ |
||
| u2 = | −dn+1(n + 1)[α β1] | |
| −dn+1(n + 1)[α β2] | ||
| f2 = | −d−2(n + 1)(2dn+1 + α2dn−1(n − 1) + |
|
| −2d(β1(β1 + 2β2) + β2(2β1 + β2)) − dn−1 n(n + 1) | ||
| −d−2(n + 1)(2dn+1 + α2dn−1(n − 1) + |
||
| λ1 = | ω | Γ1 |
| λ2 = | ω | Γ2 |
| p2 = | 0 | ∂1Ω2 |
| p2 = | 0 | ∂2Ω2 |
| p1 = | −ω | Ω1 |
| u1 = | Ω1 | |
| ∑j∈𝓢̂1 λj = | 2ω | Ω1 |
| f1 = | Ω1 |
D.2 Exact solutions for the 2D/3D validation
Analogously to the previous case, we decompose the three-dimensional matrix Ω2 into nine subdomains, i.e.,
The distance function d2(x) from Ω2 to Ω1 is now given by
and the bubble function ω(x):
In Table 8, we show the exact solutions for all the variables of interest. Once again, a value of n = 1.5 is adopted for this validation.
Exact solutions for the 2D/3D validation.
| p2 = | dn+1 + ω d | |
|---|---|---|
| Ω2 ∖ |
||
| u2 = | −dn−1(n + 1)[α β1 γ1] | |
| −dn−1(n + 1)[α β1 0] | ||
| −dn−1(n + 1)[α β1 γ2] | ||
| −dn−1(n + 1)[α 0 γ1] | ||
| −dn−1(n + 1)[α 0 γ2] | ||
| −dn−1(n + 1)[α β2 γ1] | ||
| −dn−1(n + 1)[α β2 0] | ||
| −dn−1(n + 1)[α β2 γ2] | ||
| f2 = | ||
| −α−2ω dn+1(n + 1)2 − α−2ω dn+1(n + 1) | ||
| λ1 = | ω | Γ1 |
| λ2 = | ω | Γ2 |
| p2 = | 0 | ∂1Ω2 |
| p2 = | 0 | ∂2Ω2 |
| p1 = | −ω | Ω1 |
| u1 = | Ω1 | |
| ∑j∈𝓢̂1 λj = | 2ω | Ω1 |
| f1 = | Ω1 |
© 2023 Jhabriel Varela, Elyes Ahmed, Eirik Keilegavlen, Jan M. Nordbotten, and Florin A. Radu, published by De Gruyter
This work is licensed under the Creative Commons Attribution 4.0 International License.
