1. Introduction
Oil production or fluids injection induce stress variations within the reservoir, driving compaction and surface’s subsidence with severe effects in well-casings, cap-rock, fault reactivation, as well as adverse environmental impact. We then require reliable computations to estimate these variations and their influence. Towards that end, we must pair flow simulations with poroelasticity that substantially rises memory consumption and processing time, e.g., the mechanics component can take as much as 90% of the total runtime [
1,
2,
3,
4,
5,
6,
7,
8]. Domain decomposition procedures have been intensively analyzed to overwhelm the central processing unit (CPU) load related to these predictions, and they involve splitting the whole domain into smaller subdomains in such a fashion that each part can be tackled individually to use parallel computing.
Poroelasticity is the underlying theory to estimate compaction and related risks, i.e., surface subsidence and near-wellbore issues induced by hydrocarbons production [
7,
9,
10,
11]. Rice and Cleary [
12] reformulated Biot’s theory [
13] to express the equations regarding coefficients that are more suitable for practical applications [
14,
15]. The finite element method (FEM) is the most popular tool to solve the poroelasticity equations. Indeed, Lewis and Schrefler [
16] used such techniques to 1- and 2-D consolidation problems. Phillips and Wheeler [
17,
18] showed the theoretical convergence of 2-D models that pair both continuous (CG) and discontinuous Galerkin schemes for the displacements with mixed finite element (FE) spaces for flow. Sanhu and Wilson [
19] published an early FE application to 3-D consolidation that was reviewed by Yokoo [
20]. Ghaboussi and Wilson [
21] considered fluid’s compressibility as well as the stability of CG schemes [
22]. Zienkiewicz et al. [
23] incorporated solid grains’ compressibility while Booker and Small [
24] performed the stability analysis in time. Zienkiewicz and Shiomi [
25] tackled various CG approaches for subsurface consolidation. Mixed finite elements [
17,
18], reduced quadrature rules [
26], and penalty factors [
27], which were analyzed and applied for tackling incompressibility in elasticity, was purposed for compaction computations coupled with the incompressible fluid flow. Many CG schemes have been used for various practical poromechanics applications. The reader can find particular models in petroleum engineering in [
6,
28,
29,
30].
MFEM is a reliable means to enforce a weaker continuity requirement at the boundary of subdomains in which we use various grids, e.g., non-matching. The technique permits pairing elasticity and poroelasticity on multiple domains, which is a need for flow combined with geomechanics. The reader can find a brief intro to MFEM, its evolution, and traditional references in [
1,
2,
3,
31]. Bernardi et al. proposed MFEM for Poisson’s equation to express a relaxed continuity requirement across subdomain boundaries that we use either non-conforming grids or distinct variational forms [
32]. We summarized historical and recent developments regarding MFEM in Florez [
1]. For instance, we refer the reader to papers in [
33,
34,
35,
36,
37] for applications for elasticity with curved surfaces in solid mechanics, and dual spaces [
38,
39]. Also, Florez [
1,
40] introduced a mortar scheme able to glue subdomains on curved interfaces geometrically modeled using NURBS entities, i.e., curves and surfaces. His paper included computer implementation details.
We introduced in [
41] a reservoir geometry rebuilding algorithm that generates matching hexahedral grids that we can employ for the application herein. The simplicity of the earlier procedure advises remodeling the reservoir in an equal way but letting non-matching grids in the various reservoir levels as we depict in
Figure 1. We thus reviewed our method in [
40,
42] to enhance its execution and to enable distinctive grids in the overburden and underburden compared with the pay-zone level. We present herein an MFEM that sticks these three non-conforming subdomains ensuring continuous displacements over the given boundaries that we geometrically represent via NURBS entities. We determine the SPP to require continuity among subdomains, and we further decouple the SPP with the Dirichlet-Neumann domain decomposition method (DNDDM). We briefly introduce MFEM, and we also discuss on calculating the mortar projector including all algorithms in pseudo-code. The preliminary numerical results confirm that MFEM can undertake problems whose size is of practical interest yet at the whole reservoir field.
The remaining sections are:
Section 2 defines the mathematical model and finite element formulation for single-phase flow coupled with poroelasticity.
Section 3 comprises MFEM and emphasizes on computing the projection matrix, and thus
Section 4 includes uncoupling procedures applying domain decomposition approaches.
Section 5 involves numerical examples while the three last sections encompass conclusions, path-forward, and acknowledgments sequentially.
2. Governing Equations and Finite Element Formulation
We examine the equations that govern homogeneous linear isotropic poroelasticity and introduce their FE formulation. We exclude all particulars herein for the sake of conciseness. The reader can find a more comprehensive development in [
2,
3,
4,
6,
16,
22]. We consider a bordered domain
and its outskirt
, and
is a non-deteriorated, almost even conforming triangulation of
based on hexahedra. We can show that [
43], for flexible porous media, the single-phase flow model develops from the incompressible equation of continuity and Darcy’s law, this is:
where
is a particular algorithm porosity,
means the absolute permeability, fluid properties are the dynamic viscosity,
, and
, the density. Gravity acceleration constant is
g,
p refers to fluid pressure while
q expresses sources and sinks, e.g., injector and producer wells. The resulting model porosity,
appears as:
above
means Biot’s constant,
refers the displacement vector while
comprises the original volumetric strain,
M means Biot’s modulus [
15],
and
refer to initialization quantities. The common boundary conditions (BCS) for
p involve no-flow or homogenous Neumann, i.e.:
we also require the original pressure in the volume, and
corresponds to the outer unitary normal vector. For mechanics, we consider the quasi-steady equilibrium equation, i.e., we neglect velocity variations:
where the stress tensor is
,
means the body forces vector, e.g., gravity. BCS often concern given tractions on the part of the border. We can decompose boundary conditions in categories such as Dirichlet, i.e.,
, and Neumann,
, i.e., where we prescribe the remote tractions. Biot’s poroelasticity theory and Hooke’s law determine
by:
where the elastic moduli for linear isotropic elasticity is
,
means the Kronecker delta whereas
,
G, are the Lamé constants, and
denotes the fourth order identity tensor. The symmetrical gradient or strain tensor,
gives:
We develop a weak form by replacing Equation (
2) into Equation (
1), and we multiply by a test-function
. We then integrate across the domain and employ the theorem of divergence:
We obtain Equation (
4) weak form similarly by testing against the virtual displacement,
, which yields:
where
means tractions exerted in the boundary. The above equation resembles the virtual work theorem. We often take the continuous FE space as a finite-dimensional subset of the Sobolev spaces [
44], thus:
where
refers to the space of polynomials of total degree less than or equal to
k,
denotes test functions that are continuous across element boundaries. We expand primary variables within the element
e, i.e.,
p and
, as nodal quantities multiplied by shape functions [
45,
46]:
where
and
are shape function matrices:
above
means given element’s number of nodes while nDOF denotes the number of degrees-of-freedom that ties the space dimension,
n. The engineering strain
corresponds to:
where the linear operator
becomes:
We finalize replacing Equation (
5) into (
8) and using (
7), which renders the FE formulation, i.e.:
the above matrices expressions are:
We can arrive at the loosely coupled method in various ways. We present an instance in Equation (
16) whereas we first solve displacements by considering pressures from the prior time step. Next, we update them using the latest displacements:
and the new matrixes are:
where the factor of implicitness is
while
represents the timestep size. We remark that an iterative coupling scheme relies on the loosely coupled by internally iterating to refresh lagged quantities [
6,
22]. We realized in [
1,
2] that flow simulations, in particular, multi-phase, demand mass-conservation. Indeed, besides we discretize herein the single-phase flow equation with a Continuous Galerkin (CG) formulation is possible to post-process to produce mass-conservative fluxes locally as indicated in [
47]. We believe that CG produces meaningful pressure fields that we can employ for coupling with the mechanics.
3. The Mortar Approach
The fundamental goal here is to advance MFEM to stick rounded interfaces such as the one that
Figure 1 displays. We can describe MFEM, i.e., for linear isotropic elasticity via the following bilinear forms,
a and
, that become [
36,
37]:
where
enforces for the pasting condition across subdomain boundaries, whereas we require the displacement jump, i.e.,
, to be zero in a “weak” or integral sense:
in Equation (
19)
is the mortar space while
and
are the weighting and Lagrange multiplier spaces, respectively. Let
be a two-parameter mapping that describes a given mortar geometry, i.e., NURBS surface, where the parameters
, correspond to orthogonal directions related with the surface in the reference or computational plane,
, a rectangle whose image becomes the surface in 3-D. Let
be a conforming triangulation of the reference plane,
, whose mapping or image depicts the mortar as a discrete geometrical object, e.g., a quadrilateral surface. We also take the mortar space as a finite-dimensional subset of the continuous set of Sobolev spaces:
where
means the set of polynomials whose total degree is less than or equal to
k,
denotes continuous test functions across
edges. We can write (
19) in matrix form as:
The above equation is the so-called saddle-point problem (SPP). Notice that subdomains are only glued by the Lagrange multiplier
, if these latter were known (i.e., for elasticity, they are the uncomputed tractions on the subdomain interfaces), then we can unlink the system (
21). For further details, please refer to the following references [
1,
48,
49].
We denote the rectangular matrixes
, “projectors” because they permit us to map onto and from the mortar space. For 3-D problems, they are defined by:
where
represents the surface’s normal vector metric and
is the vector of parameters. We already presented implementation details in Florez [
1]. Herein we extend the previous treatment, and we keep some of the algorithms already reported for the sake of completeness. MFEM requires a pre-processing step to ensuring that most computations are explicit [
1], which implies possessing certain data structures. Such a step encompasses the geometrical modeling of the mortar interfaces, and surface and volume mesh generation. Algorithm 1 depicts high-level steps that preprocessor executes, where
and
are the numbers of mortar entities and subdomains, and the mortar and non-mortar quadrilateral grids, in the physical plane, are represented as
and
separately. The function
yields the argument’s number of elements. We denote non-mortar and mortar sides as NM- and M-side respectively. For instance,
Figure 2 depicts the topology that we are interested in, i.e., three non-conforming levels. We anticipate the level in the middle, i.e., the pay-zone, to be the finer one. The figure also exhibits a sample mortar surface that is an interpolation NURBS bicubic surface. The color map resembles the
surface’s parameter.
In Algorithm 1 mesh generation begins with surface discretization. Over mortar interfaces,
, though, we have to save both physical
and parametric
coordinates. Having the former, i.e.,
, bypasses determining them [
1]. Before producing a volume mesh on the subdomains, we enforce preservation of the NM-side’s surface grid points. We then map the surface-to-volume numeration. We will require such mapping table while we are assembling the global SPP or passing messages for the decoupled case. The table designates the global degrees-of-freedom (DOF) that we have to glue (notice
Figure 3).
Algorithm 1: Preprocessing high-level steps. |
|
Algorithm 2 outlines the high-level steps to compute the projector while
Figure 4 reveals the mortar and non-mortar grids in the computational plane. These grids correspond to the problem that
Figure 2 describes. To obtain a given point’s (PT) hat function (HFNT), we first compute the compact-support (CS),
, of that PT, i.e., the array of all elements (ELMT) that include the given PT. For instance, for quadrilateral meshes, the number of ELMT in the CS can be one, two or four. ELMT in the corners, in the boundary and interior ones (excluding corners), explain these numbers.
Figure 4 also depicts the HFNT, that exists only on the CS of the given PT. Notice that the CS has an associated convex-hull, outside of this box the HFNT is nil. Steps 2 and 6, in Algorithm 2, correspond to computing these CS and their HFNT.
Algorithm 2: High-level steps to compute the projector. |
|
Now, the conditional in step 9 intersects the convex hulls of the given CS. Notice that the intersection can be empty.
Figure 5 represents the intersection in which the product of HFNT exists. The volume beneath this surface is the given projector entry, according to Equation (
22). To accurately compute the integral, we exploit the resulting box in four quadrilateral elements as shown, and we then use a six order Gauss-Legendre quadrature rule to evaluate every sub-integral. This latter provides enough accuracy for our purposes. Since the computational burden to calculate the projector lies in crossing CS, then we can implement an efficient range-search algorithm based on quadtrees [
1,
40] to expedite this task.
Finally,
Figure 6 depicts the mortar projector as a matrix, for the problem described in
Figure 5, its dimensions are (11 × 11) × (9 × 9). It is indeed, in this case, a symmetric matrix. The highest values correspond to those entries where the mortar and non-mortar points are at their closest locations. The mortar projector resembles a
projector which permits mapping onto the mortar space by using the matrix
. Also, mapping to a non-mortar side relies on another method that we particularly described in [
1].
4. Domain Decomposition Techniques
Domain Decomposition Methods (DDM) are robust algorithms to determine sizable problems on high-performance computing machines. These techniques partition the global domain into various subdomains to determine the overall answer through iterative-coupling, i.e., solving the unknown BCS on the interfaces of them [
5,
50]. There is vast literature relating these schemes, but we present a simple foundation herein for the sake of completeness. Bjørstad and Widlund [
51], Bramble et al. and Marini and Quarteroni [
52] examined the Dirichlet-Neumann scheme (DNDDM), for instance, among other scholars [
53].
Let
be a symbolic differential linear operator, i.e., Laplace’s operator for example. The DNDDM scheme comprises determining a sequence of problems in the proper chain that demands a coloring tool (see
Figure 7). We paint the Dirichlet subdomains in white color while Neumann’s are black. Notice too that the border separating subdomains is
. After guessing initial values for the primary variables over
, i.e., we must provide
, we can then resolve the problem on the white subdomains (type Dirichlet) that corresponds to step 1 in Equation (
23).
We denote the primary variable “displacements” while “tractions” correspond to their orthogonal derivative in the border. We calculate the former on the boundary
after determining the first step on the white subdomains. We then transfer them to complete the second step on the Neumann subdomains. On them, since we know the tractions, we thus determine the displacements to update them for the next iteration. We relax both quantities to improve the convergence where the loosening parameters
and
belong to I =
as shown in Equation (
24):
We repeat scheme (
23) until the traction residual, namely,
, in the interface lies below a given tolerance. We also fix an allowable number of iterations to stop the process in the event of divergence. It happens that this approach entails at minimum a two-color palette for coloring or yet more [
5]. We lose parallelism since black subdomains need to wait for the whites to report their tractions. We might provide an initial-guess for tractions to alleviate this difficulty, but it is not simple for most situations. We can easily produce such guess for the multiplier
by determining a coarser-run that intends obtaining the answer over a coarser grid and then interpolating
on
. We refer the reader to [
54,
55] for more details.
5. Numerical Cases
The author coded these FE formulations into a parallel C++ application termed “Integrated Parallel Finite Element Analysis (IPFA)” whose principal features were described in [
2,
4,
5]. We ran all samples herein on a MacBook Pro laptop armed with an Intel(R) Quad-Core(TM) i7-7920HQ CPU @ 3.1GHz and sixteen GB of RAM. Notice that, we attain parallelism because the CPU is multi-core. Our iterative linear solver, i.e., conjugate gradient, uses incomplete LU decomposition (ILU) as the preconditioner, and we also use a postprocessor called “LogProc” [
56]. Currently, both IPFA and LogProc are proprietary applications [
56], but a free LogProc’s community version will be available for download. We are also considering to offer IPFA as an open-source project. We are still defining a timetable to achieve that, and we thus refer readers interested to the website listed in the above reference.
Our continuous shape functions for the spatial discretization are Lagrange polynomials, , as well as mortar entities . We only consider piecewise linear polynomials and the following numerical values for for all DN examples below. We include four sample problems next. The first one tackles single-phase steady-state flow where we prescribe the pressure field to allows us to compute errors in for MFEM via SPP and DN. The second one corresponds to the near-borehole section that we solve both serially, in parallel, and with MFEM. We provide time data and relative errors on the displacements to compare these DD schemes. The last two examples encompass coupled flow and mechanics computations on realistic reservoirs to prove the application of the recommended approaches to industrial size problems. We also estimate relative errors on displacements for comparison purposes.
5.1. Example 1: Manufactured Steady-State, Single-Phase Flow
The example encompasses a fabricated problem where we propose a pressure field, and then plug it into the continuity equation to obtain the reciprocal source term. The problem’s partial differential equation is as follows:
we consider a unitary cube as the domain with Dirichlet BCS. We assume the pressure to be:
knowing
p allows us computing the error of our MFEM approximation.
Figure 8 shows the three subdomains whose meshes are hexahedral.
Table 1 summarizes the metrics of these meshes. A NURBS interpolation surface represents the mortar geometry, whose finite element space corresponds to a tensor-product grid of (20 × 20) linear piecewise rectangular elements. We employed a frontal-direct solver to determine the global SPP (
21) [
1]. The results that
Figure 9 depicts have good accordance with the assumed pressure field. The
error is 5.68 × 10
in the pay-zone level for instance.
Snapshots in
Figure 10 display unfolding matching steps for DN, where the iteration levels progress in reading order as indicated. The discontinuities among subdomains are evident in the three first snaps. The DNDDM used twelve internal iterations to lessen the residual in tractions below 10
.
Figure 11 pictures the final parallel answer, which agrees well with the above-assumed pressure. The
-error for the Dirichlet subdomain is about 3.60 × 10
, and for the Neumann improves to 4.05 × 10
. The resulting error values are small related to the magnitude of
p.
We now proceed to compare the DNDDM and MFEM for a simple manufactured 2-D problem. We now assume for Equation (
25) the following pressure field:
and the domain is the unitary square with pure homogeneous Dirichlet boundary conditions.
Figure 12 depicts the meshes partitioning and coloring to solve the problem with both schemes. In the left-hand side, we employ for the DNDDM a global conforming mesh that encompasses 9600 quadrilateral elements, and we partition it into three almost even subdomains as depicted. We color the center domain as Neumann and the remainders as Dirichlet type. Similarly, for MFEM we consider three subdomains, two of them (the top and bottom ones) are triangular meshes while the one in the middle is a regular Cartesian quadrilateral mesh as shown in the right-hand-side of the figure.
Table 2 describes the type, the number of elements and points of each mesh from top to bottom. We first focus on the MFEM problem and thus the mortars as geometrical entities consist of two B-Splines interpolants that we constructed by interpolating a sinusoidal wave as revealed. We employed thirty-two quadratic mortar elements per curve to stick these three subdomains by solving the global saddle-point problem (
21) with a direct frontal solver [
57]. The results depicted on
Figure 13 reproduce the analytical solution, we thus also display the absolute error against the later, which is on the order of
. The error increases towards the interfaces as expected.
We tried to manufacture a good benchmark problem to compare these schemes, but the DN subdomains have a more refined mesh that the MFEM ones. However, for the topological point of view, the meshes are equivalent. MFEM offers the flexibility of gluing nonmatching interfaces while DN operates over conforming ones. We wish to emphasize that these schemes are rather complementary than competitive. Indeed, we foresee that DN will help us decouple MFEM’s SPP for realistic 3-D problems via the mortar mappings that we presented in [
1].
Figure 14 displays pressure snapshots that correspond to four different Dirichlet-Neumann iteration levels, namely 2, 4, 6 and 8, evolving from left-to-right and top-to-bottom. We employ the mesh in the left-hand-side of
Figure 12. We provide no initial guess for pressure, which explains the mismatch in the first snapshots. Notice then how the process to match up those subdomains could reduce discrepancies in just a few iterations. As for stopping criterion, we require that the residual in the tractions in the interface lies below a given tolerance. For this problem, we spent sixteen iterations to achieve a traction residual smaller than
.
The DN final solution that we obtained after merging the subdomain results reproduces the analytical solution as well but exhibits larger errors compared to the MFEM, see
Figure 15. Indeed, the error triplicates in particular in the bottom interface among subdomains. Part of the issue may be the fact that we are partitioning in the area where the pressure experiences strong gradients, as smarter partitioning strategy might lead to better results but with a detrimental in balance loading.
Finally,
Table 3 compares timing data for both schemes. We employ multi-threading assembling for all subdomains. MFEM duplicates DN preprocessing time since we need to compute the mortar projectors and still need communication tables. However, since we employ the SPP approach for MFEM, we do not need to match subdomains iteratively as DN does, which explains the time savings. The DN meshes were finer than MFEM ones, and thus it takes longer to both assemble and solve those subdomains problems. The next example will cover on detail the overhead associated with DN. There is a clear winner for this benchmark problem, but results are relative, we repeat that these are complementary methods and not competitive.
5.2. Example 2: 3-D Near Wellbore Problem
We revisit herein the well-known Bradley’s problem that we discussed already in [
1] (see
Figure 16). We consider homogeneous linear isotropic elasticity.
Figure 16 depicts the problem geometry and BCS for a 2-D case where we prescribed pressure in the borehole, symmetrical displacements becomes zero on the bottom and left surfaces while we enforce the in situ or remote stresses in the remaining surfaces, as depicted in the figure. The hole’s radius is 5.0 inches, and the square length is 25.0 inches.
Table 4 below shows the mechanical properties, pressure unit is MPa, and values for BCS where we divide the parameters by the mean remote stress, i.e.,
, to normalize them. We extend Bradley’s problem to 3-D as proceeds. The geometrical entity encompasses a cube at the origin with a side of 25.0 inches, and we subtract an aligned cylinder whose radius is 5.0 inches as pictured in
Figure 17. Additionally to the 2-D conditions above, we axially enforce, i.e.,
z-direction, no
displacement in the bottom face while the top surface is traction-free, as displayed in the figure.
Our objectives with this example are twofold. First, we compare the global serial solution with its parallel counterpart for a conforming mesh case. Second, we employ MFEM to glue nonmatching interfaces that discretize the problem accordingly. We compare the runtimes of the above procedures. We first tackle Bradley’s problem serially with a global conforming mesh that has 98,304 elements and 105,105 points. Assembling the stiffness matrix takes 3 min, and 14 s and solving the sparse system last 39.5 s. We use ILU as preconditioner and CG as an iterative solver which spent 283 iterations to achieve a tolerance of 9.315 × 10
.
Table 5 summarizes the overall runtime for getting a solution. The global solution renders the displacements to lie within the following intervals:
= [−0.42786, 0],
= [−0.01154, 0.103155], and
= [0, 0.188838]. We now employ the DNDDM over a pair of processors as depicted in
Figure 18, which also shows on the right-hand-side the resulting parallel solution obtained after executing ten internal iterations to attain traction residual of 4.8 × 10
. We number the conforming subdomains starting at the bottom,
Table 5 compiles timing and partitioning data.
Snaps in
Figure 19 show the DN matching process evolution, where we illustrate the
displacement. We show the first four iterations in reading order. This parallel answer matches its serial counterpart which we gauged against commercial software, and we obtained excellent accordance [
2] (see
Figure 18 for instance). Indeed, for example, the relative error in the
displacement against the global solution is less than 0.14%. For the other two displacements, namely,
and
, the errors are 0.19% and 1% respectively. We can claim that the parallel solution reproduces the global solution. However, its runtimes are disappointing.
Table 5 highlights the fact that the parallel solution took more than twice the time spent by the global one. The Neumann subdomain contributed the longer runtime. Besides we solve smaller subdomain problems, i.e., the sparse systems are about a half in rank, we must match the subdomain which takes up to 10 iterations. We solve the smaller matrixes that many times. Also, for the Neumann’s subdomain, we compute another solution by switching the boundary conditions in the interface to Dirichlet type, so that we can overrelax the tractions as indicated in Equation (
24). Bottom line, Neumann’s subdomain acts as preconditioner and thus requires two local solves per iteration. This overhead explains the additional runtime. We must point out that partitioning and iterative coupling by the BCS makes plausible to solve huge problems that we cannot tackle in a single computer, though.
We undertake Bradley’s 3-D problem with non-conforming grids as follows. Hence, we radially divided the domain into two subdomains as shown in
Figure 20. The inner subdomain has a finer mesh than the outer one. Both are tensor-product hexahedral meshes.
Table 6 summarizes information concerning these subdomain. In the table, the mesh’s size refers to the tangential times radial times the axial number of points respectively. We consider two successively refined meshes for MFEM purposes via SPP.
Figure 20 depicts the case with the finer meshes. The table also includes the number of elements and points on those subdomains as well as the runtime to assemble the local stiffness matrixes.
We geometrically describe the mortar plane by using a NURBS circular surface as pictured in
Figure 20. The pseudo-color map presents the
variable of the entity [
2,
41]. The mortar grid resembles a regular Cartesian mesh of size
. We consider piecewise linear rectangular elements.
Table 7 presents relevant information concerning these MFEM computations. We include details such as the preprocessing time, i.e., computing the projectors and preparing BCS, the time that the direct solver employed to solve the SPP, and also the overall runtime for MFEM. We observe that Case 1 attains proportional
-error to the DN case but with a substantial speedup. If we refine the subdomain further, as Case 2, we reduce the error in half, i.e., commensurable to DN, with a coarser mesh in the outer ring, but runtime is similar to the global case. Indeed, the SPP becomes highly ill-conditioned, and the direct frontal solver [
57] struggles to solve it. However, the fact of gluing these nonconforming interfaces is what makes MFEM attractive. It is also possible to improve the direct solver performance by using dual mortar spaces that yield a diagonal projector or to decouple the SPP via DD-schemes. Wohlmuth [
38] insists that her procedure is superior because the pairing condition is simpler to accomplish because fewer elements are in the compact-support of the basis functions. In such a scheme, we can represent the mortar projector by a diagonal matrix [
37]. Finally,
Figure 21 depicts these MFEM solutions that look relatively smooth.
5.3. Example 3: Synthetic 3-D Reservoir
This example encompasses a realistic reservoir geometry including a pay-zone of size (20 × 20 × 10) and its surroundings as depicted in
Figure 22. We enforce far-field BCS to all cuboid faces, i.e., no displacement in the orthogonal direction to the given face, except in the top one, which is traction-free. For the sake of simplicity, we assume a constant pressure drop of 1000.0 Psi in the pay-zone, which resembles a reservoir compaction and subsidence scenario [
58,
59]. The global conforming mesh, see
Figure 22, consists of 13,824 hexahedral elements. We consider isotropic linear elasticity with constant mechanical properties, for both the reservoir and its vicinity thus
Ksi and
.
Figure 23 depicts half of the domain as cut-away to show the resulting vertical displacement field. The deformation pattern is the standard blue compaction dome in the cap and the red build up in the base. The displacement is not symmetric, though. We employ this global conforming solution to benchmark our MFEM procedure.
We now partition the conforming case in the fashion that
Figure 24 depicts. We keep the same discretization in the pay-zone level, but we have aerially coarser overburden (12 × 12 × 8) and underburden (12 × 12 × 6). The mesh was attracted towards the pay-zone as usual. Since this is a graded mesh, we also have to grade the mortar space mesh, in the parametric plane, in a similar manner as
Figure 25 depicts.
Table 8 compiles information concerning these subdomains. We geometrically represent the mortar spaces using NURBS surfaces as in the earlier example.
The mortar space considers a graded tensor-product grid of 25 × 25 linear piecewise rectangular elements.
Figure 26 presents the MFEM answer via SPP. This solution is a suitable approximation to the conforming one, notice that the maximum compaction is within 6% while the build-up has a 7% relative error. These results are striking besides the mesh in the overburden and underburden have four times fewer elements than their equivalents in the conforming case.
Alternatively, successive snapshots, in reading order, in
Figure 27 show the evolution of the DN pairing workflow,
is the displacement depicted; this decoupled solution agrees very well with both the serial and the one based on SPP.
Figure 27 depicts the answer achieved after completing 12 iterations to attain traction residual less than 10
. This last snapshot looks relatively smooth.
5.4. Example 4: Realistic Synthetic Reservoir
We now modify the third example to tackle a more realistic reservoir geometry as shown in
Figure 28. We employ the loosely coupled formulation (
16). The isotropic porosity and permeability fields arise from the SPE10 related project (see
Figure 29). The fluid viscosity equals 0.0133 cp while the entire compressibility is 1.4 × 10
Psi
. We consider no gravity for both flow and elasticity equations. We retain the same mesh sizes, mechanical properties, and BCS as the example before.
We consider two vertical producer wells in the corners, with prescribed bottom hole pressure, as depicted in
Figure 30. The initial condition comprises a uniform pressure of 1000.0 Psi in the entire pay-zone while we assume
p in the wells as 0.0 Psi that resembles a depletion situation. We also consider no-flow BCS on all reservoir surfaces. The initial displacement field is zero as well.
Figure 30 depicts, in reading order, pressure evolution snapshots after 0, 20, and 40 years of production. We analyze these last two snapshots for MFEM purposes.
Figure 31 depicts pictures with the unfolding of the vertical displacement, for the conforming mesh. A compaction dome develops just above the section where the most significant pressure drop occurs. The deformation-pattern is again the usual one where a compaction skull sits on the head while a build-up happens in the base of the reservoir.
We repeat the same numerical experiment, this time by solving the SPP, and we observe that the resulting MFEM solutions once again accord very well with the conforming case as depicted in
Figure 32. The discrepancies between the conforming and MFEM solutions are within 6%.
6. Concluding Remarks
We introduced herein a comparably comprehensive MFEM scheme on curved interfaces described using NURBS surfaces. We explained in detail the steps that are necessary to compute the mortar projector, and we also included their algorithms in pseudo-code. The examples considered a steady-state single-phase flow problem, Bradley’s near-borehole problem, and two coupled flow and mechanics simulations with a realistic reservoir geometry. These last two field level compaction computations showed that the suggested method has the potential to tackle problems of practical interest.
We benchmarked these schemes for a topologically equivalent 2-D problem and we find out about the pros and cons of every approach. MFEM offers the flexibility of gluing nonmatching interfaces while DN operates over conforming ones. DN must match the subdomains in an iterative fashion which consumes CPU. Regarding runtime, MFEM certainly speeds up a given problem where we can glue coarser meshes in those areas with smaller gradients. On the other hand, MFEM via SPP requires solving an ill-condition problem that may not be tractable in challenging 3-D problems. We thus emphasize that these schemes are rather complementary than competitive.
We demonstrated that the DNDDM could decouple the global SPP that is imperative to continue this methodology to 3-D domains with a substantial number of DOF, where solving the latter is not feasible at all. We found that often parallel solution runtimes are disappointing when compared with the serial ones. Indeed, the Neumann subdomain contributed the longer runtime since we compute another solution by switching the BCS in the interface to Dirichlet type so that we can overrelax the tractions. Bottom line, Neumann’s subdomain acts as preconditioner and thus requires two local solves per iteration. We must point out that partitioning and iterative coupling by the BCS makes plausible to solve huge problems that we cannot tackle in a single computer, though. This study reveals new prospects to advance parallel codes for reservoir simulation linked with poroelasticity.