Robust linear domain decomposition schemes for reduced non-linear fracture flow models

TL;DR

This paper introduces two robust linear domain decomposition algorithms for simulating non-linear fracture flow in porous media, utilizing the L-scheme for linearization and multiscale flux basis for computational efficiency, with proven stability and convergence.

Contribution

The paper develops two new algorithms, MoLDD and ItLDD, that efficiently handle non-linear fracture flow problems using linearization and precomputed flux bases, advancing computational methods in this field.

Findings

Algorithms are stable and convergent under optimized parameters.

Extensive numerical tests confirm theoretical predictions.

Pre-computed flux basis reduces computational complexity.

Abstract

In this work, we consider compressible single-phase flow problems in a porous media containing a fracture. In the latter, a non-linear pressure-velocity relation is prescribed. Using a non-overlapping domain decomposition procedure, we reformulate the global problem into a non-linear interface problem. We then introduce two new algorithms that are able to efficiently handle the non-linearity and the coupling between the fracture and the matrix, both based on linearization by the so-called L-scheme. The first algorithm, named MoLDD, uses the L-scheme to resolve the non-linearity, requiring at each iteration to solve the dimensional coupling via a domain decomposition approach. The second algorithm, called ItLDD, uses a sequential approach in which the dimensional coupling is part of the linearization iterations. For both algorithms, the computations are reduced only to the fracture by pre-computing, in an offline phase, a multiscale flux basis (the linear Robin-to-Neumann co-dimensional map), that represent the flux exchange between the fracture and the matrix. We present extensive theoretical findings and in particular, the stability and the convergence of both schemes are obtained, where user given parameters are optimized to minimise the number of iterations. Examples on two important fracture models are computed with the library PorePy and agree with the developed theory.

Figures

a](#S8.F4.sf1), we show the results for the ItLDD solver on a set of realizations of $(L_{\gamma,u},L_{\gamma,p})$. The results do not differ greatly comparing to the case of Forchheimer’s flow model. The convexity of the surface plots in all time steps is clear giving away an optimal combination of values for $L_{\gamma,u}$ and $L_{\gamma,p}$. For example, we can find minimum of $5$ iterations for $L_{\gamma,u}$ between approximately $0.73$ and $1.57$, and $L_{\gamma,p}$ between $10^{2.8}$ and $10^{3}$. Note that the parameters prescribed by the theoretical results, $L_{\gamma,u}=0.5$ and $L_

Tables

Table 1. Table 1: Results for the example of Subsection 8.1 . Top two tables correspond to solving with the MoLDD scheme, while bottom two correspond to solving with the ItLDD scheme. On the left we report the number of iterations by varying the mesh size h ℎ h for a fixed time step τ = 2 − 4 𝜏 superscript 2 4 \tau=2^{-4} , while on the right depending on the time step size τ 𝜏 \tau for a fixed mesh size h = 2 − 5 ℎ superscript 2 5 h=2^{-5} .

$h\backslash n$	1	2	3	4	5
$2^{-1}$	17	12	11	10	9
$2^{-3}$	17	11	10	9	8
$2^{-5}$	17	11	10	9	8

Table 2. Table 2: Results for the example of Subsection 8.2 reporting the number of iterations by varying the parameter β 𝛽 \beta (top) and by varying α γ subscript 𝛼 𝛾 \alpha_{\gamma} (bottom). Left tables correspond to solving with the MoLDD solver, while the right ones correspond to solving with the ItLDD solver.

$\beta \backslash n$	1	2	3	4	5
0.1	17	9	8	7	7
1	17	9	8	7	7
100	9	8	7	6	5

Table 3. Table 3: Results for the example of Subsection 8.3 reporting the number of iterations for conforming and non-conforming (coarse scale) grids on the fracture. Left table corresponds to solving with the MoLDD scheme, while the right one corresponds to solving with the ItLDD scheme.

$♯\text{cells}\backslash n$	1	2	3	4	5
8	3	3	3	3	2
16	3	3	3	3	2
64	3	3	3	3	2

Table 4. Table 4: Results for the example of Subsection 8.4 . On the left the number of iterations by varying the values of ω 𝜔 \omega . In the center when c ω subscript 𝑐 𝜔 {{c_{\omega}}} changes, while on the right for different values of r 𝑟 r .

$\omega \backslash n$	1	2	3	4	5
0.1	15	11	10	9	8
1	10	9	8	8	8
2.5	30	19	16	14	12

Equations

$${\mathbf{K}}_i^{-1}{\mathbf{u}}_i+\nabla p_i=0$$

$${\partial }_t{p}_i+\nabla \cdot {\mathbf{u}}_i=f_i$$

$$p_i=0$$

$$p_i(\cdot ,0)=p_i^0$$

$$\xi ({\mathbf{u}}_{\gamma })+{\mathbf{K}}_{\gamma }^{-1}{\mathbf{u}}_{\gamma }+{\nabla }_{\tau }{p}_{\gamma }=0$$

$${\partial }_t{p}_{\gamma }+{\nabla }_{\tau }\cdot {\mathbf{u}}_{\gamma }=f_{\gamma }+\left({\mathbf{u}}_1\cdot {\mathbf{n}}_1+{\mathbf{u}}_2\cdot {\mathbf{n}}_2\right)$$

$$p_{\gamma }=0$$

$$p_{\gamma }(\cdot ,0)=p_{\gamma }^0$$

$$-{\mathbf{u}}_i\cdot {\mathbf{n}}_i+{\alpha }_{\gamma }{p}_i$$

$${\mathbf{V}}_i:=\left\{\mathbf{v}\in \mathbf{H}(\text{div},{\Omega }_i):\mathbf{v}\cdot {\mathbf{n}}_i\in L^2(\gamma )\right\}\text{and}{M}_i:=L^2({\Omega }_i),$$

$$a_i(\mathbf{u},\mathbf{v}):=({\mathbf{K}}^{-1}\mathbf{u},\mathbf{v}{)}_{{\Omega }_i}+{\alpha }_{\gamma }^{-1}(\mathbf{u}\cdot {\mathbf{n}}_i,\mathbf{v}\cdot {\mathbf{n}}_i{)}_{\gamma },b_i(\mathbf{u},q):=(\nabla \cdot \mathbf{u},q{)}_{{\Omega }_i},c_i(p,q):=(p,q{)}_{{\Omega }_i}.$$

$$a_{\gamma }(\mathbf{u},\mathbf{v}):=({\mathbf{K}}_{\gamma }^{-1}\mathbf{u},\mathbf{v}{)}_{\gamma },b_{\gamma }(\mathbf{u},\mu ):=({\nabla }_{\tau }\cdot \mathbf{u},\mu {)}_{\gamma },c_{\gamma }(\lambda ,\mu ):=(\lambda ,\mu {)}_{\gamma }.$$

$$a_i(\mathbf{u},\mathbf{v})-b_i(\mathbf{v},p)+(p_{\gamma },\mathbf{v}\cdot {\mathbf{n}}_i{)}_{\gamma }=0$$

$$c_i({\partial }_{t}p,q)+b_i(\mathbf{u},q)=(f,q{)}_{{\Omega }_i}$$

$$(\xi ({\mathbf{u}}_{\gamma }),\mathbf{v}{)}_{\gamma }+a_{\gamma }({\mathbf{u}}_{\gamma },\mathbf{v})-b_{\gamma }(\mathbf{v},p_{\gamma })=0$$

$$\displaystyle c_{\gamma}(\partial_{t}p_{\gamma},\mu)+b_{\gamma}({\mathbf{u}}_{\gamma},\mu)-(\left\llbracket{\mathbf{u}}\cdot{\mathbf{n}}\right\rrbracket,\mu)_{\gamma}=(f_{\gamma},\mu)_{\gamma}\quad$$

$$\nabla \cdot {\mathbf{V}}_h=M_h,\text{ and }{\nabla }_{\tau }\cdot {\mathbf{V}}_{h,\gamma }=M_{h,\gamma },$$

$$(\nabla \cdot (\mathbf{u}-{\tilde{\Pi }}_i\mathbf{u}),q{)}_{{\Omega }_i}=0$$

$$((\mathbf{u}-{\tilde{\Pi }}_i\mathbf{u})\cdot {\mathbf{n}}_i,\mathbf{v}\cdot {\mathbf{n}}_i{)}_{\partial {\Omega }_i}=0$$

$$||{\tilde{\Pi }}_i\mathbf{u}|{|}_{{\Omega }_i}\lesssim ||\mathbf{u}|{|}_{ϵ,{\Omega }_i}+||\nabla \cdot \mathbf{u}|{|}_{{\Omega }_i}.$$

$$||\lambda |{|}_{\gamma }\lesssim ||{\mathcal{Q}}_{h,1}\lambda |{|}_{\gamma }+||{\mathcal{Q}}_{h,2}\lambda |{|}_{\gamma },$$

$$a_i({\mathbf{u}}_h^n,\mathbf{v})-b_i(\mathbf{v},p_h^n)=-(p_{h,\gamma }^n,\mathbf{v}\cdot {\mathbf{n}}_i{)}_{\gamma }$$

$$c_i(p_h^n-p_h^{n-1},q)+{\tau }^n{b}_i({\mathbf{u}}_h^n,q)={\tau }^n(f^n,\mu {)}_{{\Omega }_i}$$

$$(\xi ({\mathbf{u}}_{h,\gamma }^n),\mathbf{v})+a_{\gamma }({\mathbf{u}}_{h,\gamma }^n,\mathbf{v})-b_{\gamma }(\mathbf{v},p_{h,\gamma }^n)=0$$

$$\displaystyle c_{\gamma}(p^{n}_{h,\gamma}-p^{n-1}_{h,\gamma},\mu)+\tau^{n}b_{\gamma}({\mathbf{u}}^{n}_{h,\gamma},\mu)-\tau^{n}(\left\llbracket{\mathbf{u}}_{h}^{n}\cdot{\mathbf{n}}\right\rrbracket,\mu)_{\gamma}=\tau^{n}(f_{\gamma}^{n},\mu)_{\gamma}$$

$$p_{h,i}^n=p_{h,i}^{\ast }({\lambda }_{h,\gamma }^n)+{\bar{p}}_{h,i}^n\text{and}{\mathbf{u}}_{h,i}^n={\mathbf{u}}_{h,i}^{\ast }({\lambda }_{h,\gamma }^n)+{\bar{\mathbf{u}}}_{h,i}^n,1\le n\le N,$$

$$a_i({\mathbf{u}}_{h,i}^{\ast }({\lambda }_{h,\gamma }^n),\mathbf{v})-b_i(\mathbf{v},p_{h,i}^{\ast }({\lambda }_{h,\gamma }^n))=-({\lambda }_{h,\gamma }^n,\mathbf{v}\cdot {\mathbf{n}}_i{)}_{\gamma }$$

$$c_i(p_{h,i}^{\ast }({\lambda }_{h,\gamma }^n),q)+{\tau }^n{b}_i({\mathbf{u}}_{h,i}^{\ast }({\lambda }_{h,\gamma }^n),q)=0$$

$$a_i({\bar{\mathbf{u}}}_i^n,\mathbf{v})-b_i(\mathbf{v},{\bar{p}}_{h,i}^n)=0$$

$$c_i({\bar{p}}_{h,i}^n-p_{h,i}^{n-1},q)+{\tau }^n{b}_i({\bar{\mathbf{u}}}_{h,i}^n,q)={\tau }^n(f^n,\mu {)}_{{\Omega }_i}$$

$$({\bar{p}}_{h,i}^0,\mu {)}_{{\Omega }_i}=(p_{h,i}^0,\mu {)}_{{\Omega }_i}$$

$$s_{\gamma ,i}({\lambda }_{h,\gamma }^n,\mu )$$

$$s_{\gamma }({\lambda }_{h,\gamma }^n,\mu )$$

$$g_{\gamma }^n(\mu )$$

$$(\xi ({\mathbf{u}}_{h,\gamma }^n),\mathbf{v}{)}_{\gamma }+a_{\gamma }({\mathbf{u}}_{h,\gamma }^n,\mathbf{v})-b_{\gamma }(\mathbf{v},{\lambda }_{h,\gamma }^n)=0$$

$$c_{\gamma }({\lambda }_{h,\gamma }^n-{\lambda }_{h,\gamma }^{n-1},\mu )+{\tau }^n{b}_{\gamma }({\mathbf{u}}_{h,\gamma }^n,\mu )+{\tau }^n{s}_{\gamma }({\lambda }_{h,\gamma }^n,\mu )={\tau }^n(f_{\gamma }^n+g_{\gamma }^n,\mu {)}_{\gamma }$$

$$(\xi ({\mathbf{u}}_{h,\gamma }^{n,k-1})+L_{\gamma }({\mathbf{u}}_{h,\gamma }^{n,k}-{\mathbf{u}}_{h,\gamma }^{n,k-1}),\mathbf{v}{)}_{\gamma }+a_{\gamma }({\mathbf{u}}_{h,\gamma }^{n,k},\mathbf{v})-b_{\gamma }(\mathbf{v},{\lambda }_{h,\gamma }^{n,k})=0,$$

$$c_{\gamma }({\lambda }_{h,\gamma }^{n,k}-{\lambda }_{h,\gamma }^{n-1},\mu )+\tau b_{\gamma }({\mathbf{u}}_{h,\gamma }^{n,k},\mu )+\tau s_{\gamma }({\lambda }_{h,\gamma }^{n,k},\mu )=\tau (f_{\gamma }^n+g_{\gamma }^n,\mu {)}_{\gamma }.$$

$${\mathcal{A}}_{\text{DD}}\left[\begin{array}{c}{\mathbf{u}}_{h,\gamma }^{n,k}\\ {\lambda }_{h,\gamma }^{n,k}\end{array}\right]:=\left[\begin{array}{cc}{\mathbf{A}}_{L,\gamma }& {\mathbf{B}}_{\gamma }^{\text{T}}\\ {\mathbf{B}}_{\gamma }& {\mathcal{S}}_{\gamma }^{\text{RtN}}+\mathbf{I}/\tau \end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{h,\gamma }^{n,k}\\ {\lambda }_{h,\gamma }^{n,k}\end{array}\right]=\left[\begin{array}{c}{L}_{\gamma }{\mathbf{u}}_{h,\gamma }^{n,k-1}-\xi ({\mathbf{u}}_{h,\gamma }^{n,k-1})\\ g_{\gamma }^n+f_{\gamma }^n+{\lambda }_{h,\gamma }^{n-1}/\tau \end{array}\right]:={\mathcal{F}}_{\gamma },$$

$${\mathcal{S}}_{\gamma }^{\text{RtN}}({\lambda }_{h,\gamma })=\sum _{i\in \{1,2\}}{\mathcal{S}}_{\gamma ,i}^{\text{RtN}}({\lambda }_{h,\gamma }).$$

Full text

Key words: Porous medium; reduced fracture models; generalized Forchheimer’s laws; mortar mixed finite element; multiscale flux basis; non-linear interface problem; non-overlapping domain decomposition; L-scheme.

1 Introduction

Fractures are ubiquitous in porous media and strongly affect the flow and transport. Several energy and environmental applications including carbon sequestration, geothermal energy, and ground-water contamination involve flow and transport problems in a porous medium containing fractures. Typically, fractures are thin and long formations that correspond to a fast pathway along which the medium properties, such as permeability or porosity, differ from the adjacent formations (the rocks) [6, 36, 41, 19]. Specifically, the permeability of the fracture can be significantly higher than that of the host rock. As a consequence, while flow in the host rock can be well represented by the linear Darcy’s law, flow in the fractures can potentially exhibit non-linear effects. Models for such flow will thus be non-linear, but with the non-linear effects confined to specific parts of the domain, moreover, these regions are characterized by an extreme aspect ratio.

In this paper, we consider models with fracture flow represented by an unstready Forchheimer’s law  [35, 36], as an extension of the model in [2, 36]. Our approach can straightforwardly be broadened to cover viscosity models for generalized Newtonian fluids [24, 23]. We also refer to [37, 1, 49, 4, 38] for extensions to other flow models. To limit complications relating to mesh construction and computational cost, we represent the fracture as a lower-dimensional object embedded in the full domain, as introduced in [6, 14], and we refer to the resulting model as mixed-dimensional. We then apply a domain decomposation (DD) approach, with the fracture forming an interface between subdomains. The domain decomposition approach is beneficial for both modeling, discretization and the formulation non-linear and linear solvers

Considerable research efforts have been conducted to mixed-dimensional fracture models. Several numerical schemes for steady-state models have been proposed, such as the cell-centered finite volume scheme [31], the extended finite element method [21], the mimetic finite difference [8], the block-centered finite difference method [40] and the mixed finite element (MFE) methods [41, 20, 25, 15]. Herein, we discretize the generalized mixed-dimensional Forchheimer problem with a mortar mixed finite element method (MMFEM) [29, 50, 10], combined with backword Euler in time. See also [13] for a review on fracture models and discretization approaches.

While the dimension reduction reduces the number of cells necessary to represent the fracture, the computational cost in solving the discrete non-linear problem can still be significant, in particular for complex fracture geometries. This calls for the construction of efficient solvers, and domain decomposition facilitates the exploitation of both the geometric structure of the problem, and the spatial localization of non-linearities. See, for example, the application of DD to reduced Darcy [6, 33] and Darcy-Forchheimer [26] fracture models. In this work, we develop DD schemes [2, 6, 43] to solve the non-linear problem resulting from our discretization. To exploit the geometric structure of the problem, we reformulate it as an interface problem by eliminating the subdomain variables, obtaining a non-linear system to solve at each time step. That is, the resulting problem posed only on the fracture is a superposition of a non-linear local flow operator within the fracture and a linear non-local one (Robin-to-Neumann type) handling the flux contributions from the subdomains. For this problem, two schemes are proposed, both based on the so-called $L$-scheme method, a robust quasi-Newton method with a parameter $L>0$ mimicking the Jacobian  [44, 39]. Our two approaches differ in the way they handle the non-linearity, and in the degree of coupling between the fracture interface and the surrounding subdomains.

The first algorithm named the Monolithic LDD-scheme (MoLDD) employs the $L$-scheme as a linearization procedure. At each $L$-scheme iteration, an inner algorithm is used to solve the linear interface problem [27]. It can be a direct or an iterative method (e.g. a Krylov method). The action of the interface operator requires solving subdomain problems with Robin boundary conditions on the fracture. This algorithm is Jacobian-free, solving subdomain problems can be done in parallel and is later shown to be unconditionally stable. We also obtain the condition number estimates of the inner DD system, the contraction estimates and rates of convergence for the outer scheme. However, there is still a computational overhead associated with its non-local part [2, 27, 1], that is, the subdomain solves. Increasing the non-linearity strength, the number of subdomains and refining the grids all lead to an increase in the number of iterations and the number of subdomain solves.

More recently, the $L$-scheme has gained attention as an efficient solver to treat simultaneously non-linear and coupling effects in complex problems [45, 16]. Building on this idea, we propose the second algorithm, referred as the Iterative LDD-scheme (ItLDD). In ItLDD, the $L$-scheme is now synchronizing linearization and domain decomposition through one-loop algorithm [5, 17]. At each iteration it has the cost of the sequential approach, yet it converges to the fully monolithic approach. This way we reduce the computational cost as no inner DD solver is required and only a modest number of subdomain solves that can be done in parallel are needed at each iteration. This algorithm increases local to non-local cooperation and saves computational time if one process is dominating the whole problem. This approach differs from the one commonly used in DD methods for non-linear interface problems [3, 12].

The second contribution of this paper concerns the robust and efficient implementation of the two LDD schemes. The dominant computational costs in these schemes comes from the subdomain solves and, to reduce this, we use the multiscale flux basis framework from [29]. The fact that the non-linearity in the system appears within the local operator on the fracture motivates that the linear non-local contribution from the subdomains can be expressed as a superposition of multiscale basis functions [29, 28, 2], in spirit of reduced basis [48, 32, 11]. This multiscale flux basis consists of the flux (velocity trace) response from each fracture pressure degrees of freedom. They are computed by solving a fixed number of steady Robin subdomain problems, equal to the number of fracture pressure degrees of freedom per subdomain. An inexpensive linear combination of the multiscale flux basis functions then replaces the subdomain solves in any inner/outer iteration of the algorithms. This step of freezing the contributions from the rock matrices can be cheaply evaluated and easily implemented in the algorithms. It permits reusing the same basis functions to compare MoLDD with ItLDD, to simulate various linear and non-linear models for flow in the fracture and to vary the fracture permeability. In case of a fixed time step, the multiscale flux basis is constructed only once in the offline phase. Numerical results are computed with the library PorePy [34].

In Section 2 the model problem is presented. The approximation of problem (2.6) using the MMFEM in space and a backward Euler scheme in time is given in Section 3. Also, the reduction of this mixed-dimensional scheme into a non-linear interface one is introduced. The LDD-schemes are formulated in Section 4. In Section 5 and Section 6, the analysis of the schemes is presented. Section 7 describes the implementation based on the multiscale flux basis framework. Finally, we show the performance of our methods on several numerical examples in Section 8 and draw the conclusions in Section 9.

2 Model problem

Let $\Omega$ be a bounded domain in ${\amsmathbb R}^{d}$, $d\in\{2,3\}$, with Lipschitz continuous boundary $\Gamma\vcentcolon=\partial\Omega$. Furthermore, let $T$ be the final time simulation and $I\vcentcolon=(0,T)$. Suppose that $\gamma\subset\Omega$ is a $(d-1)$-dimensional non self-intersecting surface of class $C^{2}$ that divides $\Omega$ into two subdomains: $\Omega=\Omega_{1}\cup\Omega_{2}\cup\gamma$, where $\gamma\vcentcolon=\partial\Omega_{1}\cap\partial\Omega_{2}$ and $\Gamma_{i}\vcentcolon=\partial\Omega_{i}\cap\partial\Omega$, $i\in\{1,2\}$. Assume the flow in $I\times\Omega_{i}$, $i\in\{1,2\}$, is given by

[TABLE]

and in $I\times\gamma$ by the following equations

[TABLE]

where the transmission conditions for $i\in\{1,2\}$ are

[TABLE]

Here, $\nabla_{\tau}$ denotes the $(d-1)$-dimensional gradient operator in the plane of $\gamma$, $\mathbf{K}_{\gamma}$ is the hydraulic conductivity tensor in the fracture, $\mathbf{K}_{i}$ is the hydraulic conductivity tensor in the subdomain $\Omega_{i}$ and ${\mathbf{n}}_{i}$ is the outward unit normal vector to $\partial\Omega_{i}$, $i\in\{1,2\}$. The function $\xi$ is a non-linear function extending the classical Forchheimer flow to more general laws. The coefficient $\alpha_{\gamma}$ is proportional to the normal component of the fracture permeability and inversely proportional to the fracture width/aperture. The functions $f_{\gamma}$ and $f_{i}$, $i\in\{1,2\}$, are source terms in the fracture and in the matrix, respectively. For simplicity, we have imposed a homogeneous Dirichlet condition on the boundary $\partial\Omega$. Finally, $p^{0}_{\gamma}$ and $p^{0}_{i}$, $i\in\{1,2\}$, are initial conditions.

The system (2.1)–(2.3) is a mixed-dimensional model for flow in fractured porous media: the equations (2.1a)-(2.1b) are Darcy’s law and mass conservation equations in the subdomain $\Omega_{i}$, while (2.2a)-(2.2b) are generalized Forchheimer’s law and mass conservation in the fracture of co-dimension one. Together these equations form a non-standard transmission problem where the fracture system sees the surrounding matrix system through the source term ${\mathbf{u}}_{1}\cdot{\mathbf{n}}_{1}+{\mathbf{u}}_{2}\cdot{\mathbf{n}}_{2}$ in (2.2b) and the matrix system communicates to the fracture through Robin interface conditions (2.3). Note that the restriction to only one fracture is made to simplify the presentation, but the model and the analysis below can easily be extended to fracture networks [2, 41].

2.1 Assumptions on the data and weak formulation

Let $D\subseteq\Omega$. For $s\geq 0$, $||\cdot||_{s,D}$ stands for the usual Sobolev norm on $H^{s}(D)$. If $s=0$, $||\cdot||_{D}$ is simply the $L^{2}$ norm and $(\cdot,\cdot)_{D}$ stands for the $L^{2}$ scalar product. We define the weak spaces in $\Omega_{i}$ for $i\in\{1,2\}$ as

[TABLE]

where we have implicitly considered the trace operator of $\mathbf{v}\cdot\mathbf{n}_{i}$. Moreover, we introduce their global versions by $\mathbf{V}\vcentcolon=\bigoplus_{i=1}^{2}\mathbf{V}_{i}\textnormal{ and }M\vcentcolon=\bigoplus_{i=1}^{2}M_{i}.$ The mixed spaces on the fracture $\gamma$, are $\mathbf{V}_{\gamma}\vcentcolon=\mathbf{H}(\textnormal{div}_{\tau},\gamma)\textnormal{ and }M_{\gamma}\vcentcolon=L^{2}(\gamma)$. For simplicity of notation, we introduce the jump \left\llbracket\cdot\right\rrbracket given by \left\llbracket{\mathbf{u}}\cdot{\mathbf{n}}\right\rrbracket\vcentcolon={\mathbf{u}}_{1}\cdot{\mathbf{n}}_{1}+{\mathbf{u}}_{2}\cdot{\mathbf{n}}_{2}. and the functions $\mathbf{K}$ and $f$ in $\Omega_{1}\cup\Omega_{2}$ such that $\mathbf{K}_{i}\vcentcolon=\mathbf{K}|_{\Omega_{i}}$, and $f_{i}\vcentcolon=f|_{\Omega_{i}}$, $i\in\{1,2\}$. We assume:

1. (A1)

$\xi:{\amsmathbb R}\rightarrow{\amsmathbb R}$ is $C^{1}$, strictly increasing and Lipschitz continuous, i.e., there exist $\xi_{m}>0$ and $L_{\xi}$ such that $\xi_{m}\leq\xi^{\prime}({\mathbf{u}})\leq L_{\xi}<\infty$. Otherwise, we ask bounded flux in (2.1)–(2.3), i.e, ${\mathbf{u}}\in\left[L^{\infty}(\Omega)\right]^{d}$, when $\xi$ is an increasing function ($\xi^{\prime}\geq 0$), and let $L_{\xi}\vcentcolon=\sup_{|{\mathbf{u}}|\leq C_{\xi}}{\xi^{\prime}({\mathbf{u}})}$, where $C_{\xi}\vcentcolon=\sup_{{\mathbf{x}}\in\overline{\Omega}}|{\mathbf{u}}({\mathbf{x}})|$. 2. (A2)

$\mathbf{K}:{\amsmathbb R}^{d}\rightarrow{\amsmathbb R}^{d}$ is assumed to be constant in time and bounded: there exist $c_{\mathbf{K}}>0$ and $C_{\mathbf{K}}$ such that $\zeta^{\textnormal{T}}\mathbf{K}^{-1}({\mathbf{x}})\zeta\geq c_{\mathbf{K}}|\zeta|^{2}$ and $|\mathbf{K}^{-1}({\mathbf{x}})\zeta|\leq C_{\mathbf{K}}|\zeta|$ for a.e. ${\mathbf{x}}\in\Omega_{1}\cup\Omega_{2}$, $\forall\zeta\in{\amsmathbb R}^{d}$. 3. (A3)

$\mathbf{K}_{\gamma}:{\amsmathbb R}^{d-1}\rightarrow{\amsmathbb R}^{d-1}$ is assumed to be constant in time and bounded; there exist $c_{\mathbf{K},\gamma}>0$ and $C_{\mathbf{K},\gamma}$ such that $\zeta^{\textnormal{T}}\mathbf{K}^{-1}_{\gamma}({\mathbf{x}})\zeta\geq c_{\mathbf{K},\gamma}|\zeta|^{2}$ and $|\mathbf{K}^{-1}_{\gamma}({\mathbf{x}})\zeta|\leq C_{\mathbf{K},\gamma}|\zeta|$ for a.e. ${\mathbf{x}}\in\gamma$, $\forall\zeta\in{\amsmathbb R}^{d-1}.$ 4. (A4)

The Robin parameter $\alpha_{\gamma}$ is a strictly positive constant: $\alpha_{\gamma}>0$. 5. (A5)

The initial conditions are such that $p_{i}^{0}\in L^{2}(\Omega_{i})$, $i\in\{1,2\}$, and $p^{0}_{\gamma}\in L^{2}(\gamma)$. The source terms are such that $f_{i}\in L^{2}(0,T;L^{2}(\Omega_{i}))$, $i\in\{1,2\}$, and $f_{\gamma}\in L^{2}(0,T;L^{2}(\gamma))$. For simplicity we further assume that $f$ and $f_{\gamma}$ are piecewise constant in time with respect to the temporal mesh introduced.

Remark 2.1** (On assumptions).**

*The Lipschitz-continuity of $\xi$ is not true when the function $\xi$ (therefore the flux) is unbounded, as it is the case for generalized Forchheimer’s law. However, for bounded flux ${\mathbf{u}}$, this can be verified. Otherwise, this assumption can be recovered by truncating the original function $\xi$. Obviously, the solution of the truncated problem will not in general solve the original one. See, for example  [42]. *

We introduce the bilinear forms $a_{i}:\mathbf{V}_{i}\times\mathbf{V}_{i}\rightarrow{\amsmathbb R}$, $b_{i}:\mathbf{V}_{i}\times M_{i}\rightarrow{\amsmathbb R}$ and $c_{i}:M_{i}\times M_{i}\rightarrow{\amsmathbb R}$, $i\in\{1,2\}$,

[TABLE]

On the fracture, we define the bilinear forms $a_{\gamma}:\mathbf{V}_{\gamma}\times\mathbf{V}_{\gamma}\rightarrow{\amsmathbb R}$, $b_{\gamma}:\mathbf{V}_{\gamma}\times M_{\gamma}\rightarrow{\amsmathbb R}$ and $c_{\gamma}:M_{\gamma}\times M_{\gamma}\rightarrow{\amsmathbb R}$,

[TABLE]

With the above notations, a weak solution of (2.1)–(2.3) is given in the following.

Definition 2.2** (Mixed-dimensional weak solution).**

Assume that (A1)–(A5) hold true. We say that $({\mathbf{u}},p)\in L^{2}(0,T;\mathbf{V})\times H^{1}(0,T;M)$ and $({\mathbf{u}}_{\gamma},p_{\gamma})\in L^{2}(0,T;\mathbf{V}_{\gamma})\times H^{1}(0,T;M_{\gamma})$ form a weak solution of (2.1)–(2.3) if it satisfies weakly the initial conditions (2.1d) and (2.2d), and for each $i\in\{1,2\}$,

[TABLE]

Remark 2.3**.**

*In this paper we assume that a weak solution of Definition 2.2 exists. For the steady-state model and the classical Forchheimer’s law, the existence and uniqueness of a weak solution was shown in [36]. That of the linear case, i.e., $\xi\vcentcolon=0$, was studied in [33]. Through the paper, we will also consider the case of continuous pressure across $\gamma$ by letting $\alpha_{\gamma}\rightarrow\infty$ in (2.3). For this, we will use Definition 2.2 for the weak formulation changing in (2.4) to be $a_{i}({\mathbf{u}},{\mathbf{v}})\vcentcolon=(\mathbf{K}^{-1}{\mathbf{u}},{\mathbf{v}})_{\Omega_{i}}$ and $\mathbf{V}_{i}\vcentcolon=\mathbf{H}(\textnormal{div},\Omega_{i})$ for $i\in\{1,2\}$. *

3 The domain decomposition formulation

As explained earlier, it is natural to solve the mixed-dimensional problem (2.6) using domain decomposition techniques, especially as these methods make it possible to take different time grids in the subdomains and in the fracture.

3.1 Discretization in space and time

We introduce in this section the partitions of $\Omega$ and $(0,T)$, basic notation, and the mortar mixed finite element discretization of the mixed-dimensional problem (2.6) .

Let $\mathcal{T}_{h,i}$ be a partition of the subdomain $\Omega_{i}$ into either $d$-dimensional simplicial or rectangular elements. Moreover, we assume that these meshes are such that $\mathcal{T}_{h}=\displaystyle{\cup_{i=1}^{2}\mathcal{T}_{h,i}}$ forms a conforming finite element mesh on $\Omega$. We also let $\mathcal{T}_{h,\gamma}$ be either a partition of the fracture $\gamma$ induced by $\mathcal{T}_{h}$ or slightly coarser. Denote $h$ as the maximal mesh size of both $\mathcal{T}_{h}$ and $\mathcal{T}_{h,\gamma}$. For an integer $N\geq 0$, let $\left(\tau^{n}\right)_{0\leq n\leq N}$ denote a sequence of positive real numbers corresponding to the discrete time steps such that $T=\sum_{n=1}^{N}\tau^{n}$. Let $t^{0}\vcentcolon=0$, and $t^{n}\vcentcolon=\sum_{j=1}^{n}\tau^{j},\ 1\leq n\leq N$, be the discrete times. Let $I^{n}\vcentcolon=(t^{n-1},t^{n}],\ 1\leq n\leq N$.

For the approximation of scalar unknowns, we introduce $M_{h}\vcentcolon=M_{h,1}\times M_{h,2}$ and $M_{h,\gamma}$, where $M_{h,i}$, $i\in\{1,2\}$, and $M_{h,\gamma}$ are the spaces of piecewise constant functions associated with $\mathcal{T}_{h,i}$, $i\in\{1,2\}$ and $\mathcal{T}_{h,\gamma}$, respectively. For the vector unknowns, we introduce $\mathbf{V}_{h}\vcentcolon=\mathbf{V}_{h,1}\times\mathbf{V}_{h,2}$ and $\mathbf{V}_{h,\gamma}$, where $\mathbf{V}_{h,i}$, $i\in\{1,2\}$ and $\mathbf{V}_{h,\gamma}$, are the lowest-order Raviart-Thomas-Nédélec finite elements spaces associated with $\mathcal{T}_{h,i}$, $i\in\{1,2\}$ and $\mathcal{T}_{h,\gamma}$, respectively. Thus, $\mathbf{V}_{h}\times M_{h}\subset\mathbf{V}\times M$ and $\mathbf{V}_{h,\gamma}\times M_{h,\gamma}\subset\mathbf{V}_{\gamma}\times M_{\gamma}$. For all of the above spaces,

[TABLE]

and there exists a projection $\tilde{\Pi}_{i}:\mathbf{H}^{1/2+\epsilon}(\Omega_{i})\cap\mathbf{V}_{i}\rightarrow\mathbf{V}_{h,i}$, $i\in\{1,2\}$, for any $\epsilon>0$, satisfying among other properties [28] that for any ${\mathbf{u}}\in\mathbf{H}^{1/2+\epsilon}(\Omega_{i})\cap\mathbf{V}_{i}$

[TABLE]

We also note that if ${\mathbf{u}}\in\mathbf{H}^{\epsilon}(\Omega_{i})\cap\mathbf{V}_{i}$, $0<\epsilon<1$, $\tilde{\Pi}_{i}{\mathbf{u}}$ is well-defined [47] and

[TABLE]

We introduce $\mathcal{Q}_{h,i}$ the $L^{2}$-projection onto $\mathbf{V}_{h,i}\cdot{\mathbf{n}}_{i}$ and denote $\mathcal{Q}_{h,i}^{\textnormal{T}}:\mathbf{V}_{h,i}\cdot{\mathbf{n}}_{i}\rightarrow M_{h,\gamma}$ as the $L^{2}$-projection from the normal velocity trace on the subdomains onto the mortar space $M_{h,\gamma}$. Thus, for all $\lambda\in M_{h,\gamma}$

[TABLE]

can be verified if the mesh on the fracture $\mathcal{T}_{h,\gamma}$ matches the one resulting from the surrounding subdomains, or if $\mathcal{T}_{h,\gamma}$ is chosen slightly coarser [15, 9]. Note that (3.1) can be satisfied by choosing any of the usual MFE pairs. The condition (3.5) can be satisfied even if the space $M_{h,\gamma}$ is not much richer than the space of normal traces on $\gamma$ of elements of $\mathbf{V}_{h}$ [28, 29]. The fully discrete scheme of the mixed-dimensional formulation (2.6) based on the MMFEM in space and the backward Euler scheme in time is defined through the following.

Definition 3.1** (The mixed-dimensional scheme).**

At each time step $n\geq 1$, assuming $(p^{n-1}_{h,\gamma},p^{n-1}_{h})$ is given, we look for $({\mathbf{u}}^{n}_{h},p^{n}_{h})\in\mathbf{V}_{h}\times M_{h}$ and $({\mathbf{u}}^{n}_{h,\gamma},p^{n}_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ such that, for $i\in\{1,2\},$

[TABLE]

3.2 Reduction into an interface problem

Following the algorithm in [2], we reduce the mixed-dimensional scheme in Definition 3.1 to a non-linear interface one on $\gamma$. For $i\in\{1,2\}$, we let

[TABLE]

where for $\lambda^{n}_{h,\gamma}\in M_{h,\gamma}$, $({\mathbf{u}}^{*}_{h,i}(\lambda^{n}_{h,\gamma}),p^{*}_{h,i}(\lambda^{n}_{h,\gamma}))\in\mathbf{V}_{h,i}\times M_{h,i}$ solves

[TABLE]

and $(\bar{{\mathbf{u}}}^{n}_{h,i},\bar{p}^{n}_{h,i})\in\mathbf{V}_{h,i}\times M_{h,i}$ solves

[TABLE]

Define the forms $s_{\gamma,i}:M_{h,\gamma}\times M_{h,\gamma}\rightarrow{\amsmathbb R}$, $i\in\{1,2\}$, $s_{\gamma}:M_{h,\gamma}\times M_{h,\gamma}\rightarrow{\amsmathbb R}$, and $g^{n}_{\gamma}:M_{h,\gamma}\rightarrow{\amsmathbb R}$ as

[TABLE]

where $\mathcal{S}^{\textnormal{RtN}}_{\gamma,i}:M_{h,\gamma}\rightarrow M_{h,\gamma}$, $1\leq i\leq 2$, and $\mathcal{S}^{\textnormal{RtN}}_{\gamma}\vcentcolon=\sum_{i=1}^{2}\mathcal{S}^{\textnormal{RtN}}_{\gamma,i}$ are Robin-to-Neumann type operators. Consequently, the operator $\mathcal{S}^{\textnormal{RtN}}_{\gamma,i}$ is linear. It is possible to verify that the non-linear mixed-dimensional scheme (3.6) is equivalent to the non-linear interface scheme.

Definition 3.2** (The reduced scheme).**

For $n\geq 1$ and $\lambda^{n-1}_{h,\gamma}$, find $({\mathbf{u}}^{n}_{h,\gamma},\lambda^{n}_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ such that

[TABLE]

4 Robust L-type Domain-Decomposition (LDD) schemes

In this section, we propose two iterative approaches based on the $L$-scheme to solve (3.11). The first approach entails an inner-outer procedure of the form $\textit{linearize}\rightarrow\textit{solve the DD}\rightarrow update$, so that the $L$-scheme is used for the outer loop and an inner solver (direct or iterative) for the inner loop. The second approach is a one-loop procedure in which the $L$-scheme acts iteratively and simultaneously on the linearization and DD. For the presentation of the algorithms, we shall denote the time step simply by $\tau$, keeping in mind it may depend on $n$.

4.1 A monolithic LDD scheme

The monolithic LDD scheme (MoLDD) for (3.11) reads:

Algorithm 4.1** (The MoLDD scheme).**

Given $n=0$, $(\lambda^{0}_{h,\gamma},p^{0}_{h})\in M_{h,\gamma}\times M_{h}$, stabilization parameter $L_{\gamma}>0$ and tolerance $\epsilon>0$,

Do**

Increase $n\vcentcolon=n+1$. 2. 2.

Choose an initial approximation ${\mathbf{u}}^{n,-1}_{h,\gamma}\in\mathbf{V}_{h,\gamma}$ of ${\mathbf{u}}^{n}_{h,\gamma}$. Set $k\vcentcolon=-1$. 3. 3.

Do**

1. (a)

Increase $k\vcentcolon=k+1$. 2. (b)

Compute $({\mathbf{u}}^{n,k}_{h,\gamma},\lambda^{n,k}_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ such that, for all $({\mathbf{v}},\mu)\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$,

[TABLE]

while* ${\|({\mathbf{u}}^{n,k}_{h,\gamma},\lambda^{n,k}_{h,\gamma})-({\mathbf{u}}^{n,k-1}_{h,\gamma},\lambda^{n,k-1}_{h,\gamma})\|_{\gamma}}\geq\epsilon\|({\mathbf{u}}^{n,k-1}_{h,\gamma},\lambda^{n,k-1}_{h,\gamma})\|_{\gamma}$.* 4. 4.

Update the subdomain solutions via (3.7).

while* $n\leq N$.

The advantages of Algorithm 4.1 are multiple: (i) the algorithm is Jacobian-free and independent of the initialization, (ii) we can reuse the existing $d$- and $(d-1)$-dimensional codes for solving linear Darcy problem, and (iii) optimal convergence rate is obtained with a stabilization amount determined through $L_{\gamma}$.

The MoLDD scheme involves the solution of a linear Darcy interface problem (4.1) at each iteration. We introduce the linear operators $\mathbf{A}_{L,\gamma}:\mathbf{V}_{h,\gamma}\rightarrow\mathbf{V}_{h,\gamma}$ and $\mathbf{B}_{\gamma}:\mathbf{V}_{h,\gamma}\rightarrow M_{h,\gamma}$, defined as $(\mathbf{A}_{L,\gamma}{\mathbf{u}},{\mathbf{v}})_{\gamma}\vcentcolon=a_{\gamma}({\mathbf{u}},{\mathbf{v}})+L_{\gamma}({\mathbf{u}},{\mathbf{v}})_{\gamma}$, $\forall{\mathbf{u}},{\mathbf{v}}\in\mathbf{V}_{h,\gamma}$, and $(\mathbf{B}_{\gamma}{\mathbf{u}},q)\vcentcolon=b_{\gamma}({\mathbf{u}},q)$, $\forall{\mathbf{v}}\in\mathbf{V}_{h,\gamma}$, $\forall q\in M_{h,\gamma}$. (4.1) becomes

[TABLE]

which can be solved using a direct or a Krylov type method, such as GMRes or MINRes. Regardless of the choice of inner method, we need to evaluate the action of the Robin-to-Neumann type operator $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$ via (3.10), representing physically the flow contributions from the subdomains by solving Robin subdomain problems (3.8). We summarize the evaluation of the interface operator by the following steps:

Algorithm 4.2** (Evaluating the action of $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$).**

Enter interface data $\lambda_{h,\gamma}$. 2. 2.

For* $i=1:2$*

1. (a)

Project mortar pressure onto subdomain boundary, i.e., $\varphi_{h,\gamma,i}=\mathcal{Q}_{h,i}(\lambda_{h,\gamma}).$ 2. (b)

Solve the subdomain problem (3.8) with Robin data $\varphi_{h,\gamma,i}$. 3. (c)

Project the resulting flux onto the space $M_{h,\gamma}$, i.e., $\mathcal{S}^{\textnormal{RtN}}_{\gamma,i}(\lambda_{h,\gamma})=-\mathcal{Q}_{h,i}^{\textnormal{T}}{\mathbf{u}}^{*}_{h,i}(\varphi_{h,\gamma,i})\cdot{\mathbf{n}}_{i}.$

EndFor** 3. 3.

Compute the flow contribution from the subdomains to $\gamma$ given by the flux jump across the fracture,

[TABLE]

The evaluation of $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$ dominates the total computational costs in Algorithm 4.1 (step 2(b) of Algorithm 4.2). The number of subdomain solves required by this method at each time step $n\geq 1$ is approximately equal to $\sum_{k=1}^{N^{n}_{\textnormal{Lin}}}N^{k}_{\textnormal{DD}}$, where $N_{\textnormal{Lin}}$ is the number of iterations of the $L$-scheme, and $N^{k}_{\textnormal{DD}}$ denotes the number of inner DD iterations. To set up the right-hand side term $f_{\gamma}^{n}$, we need to solve once in the subdomains at each time step $n\geq 1$.

4.2 An iterative LDD-scheme

An alternative to MoLDD scheme is to let the $L$-scheme act iteratively not only on the non-linearity, but also on the fracture-matrix coupling. Additional stabilization term is then required for the inter-dimensional coupling.

Algorithm 4.3** (The ItLDD scheme).**

Given n=0, $(\lambda^{0}_{h,\gamma},p^{0}_{h})\in M_{h,\gamma}\times M_{h}$, the stabilization parameters $(L_{\gamma,p},L_{\gamma,u})>0$, and the tolerance $\epsilon>0$.

Do**

Increase $n\vcentcolon=n+1$. 2. 2.

Choose an initial approximation $({\mathbf{u}}^{n,-1}_{h,\gamma},\lambda^{n,-1}_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ of $({\mathbf{u}}^{n}_{h,\gamma},\lambda^{n}_{h,\gamma})$. Set $k\vcentcolon=-1$. 3. 3.

Do**

1. (a)

Increase $k\vcentcolon=k+1$. 2. (b)

Compute $({\mathbf{u}}^{n,k}_{h,\gamma},\lambda^{n,k}_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ such that, for all $({\mathbf{v}},\mu)\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$,

[TABLE]

while* ${\|({\mathbf{u}}^{n,k}_{h,\gamma},\lambda^{n,k}_{h,\gamma})-({\mathbf{u}}^{n,k-1}_{h,\gamma},\lambda^{n,k-1}_{h,\gamma})\|_{\gamma}}\geq\epsilon{\|({\mathbf{u}}^{n,k-1}_{h,\gamma},\lambda^{n,k-1}_{h,\gamma})\|_{\gamma}}$.* 4. 4.

Update the subdomain solutions via (3.7).

while* $n\leq N$. *

The linear problem (4.3) can again be solved by a direct or iterative method, and it requires applying the operator $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$ via Algorithm 4.2, at each time step $n\geq 1$. The advantages of the Algorithm 4.3 are: (i) at each iteration $k\geq 1$, the systems in the fracture and the rock matrices cooperate sequentially in one loop (ii) optimal convergence rate is obtained with precise stabilization parameters $(L_{\gamma,p},L_{\gamma,u})$, and (iii) existing codes for $d$- and $(d-1)$-dimensional Darcy problems can be cheaply reused.

5 Analysis of MoLDD-scheme

The complete analysis of Algorithm 4.1 will be carried out in two steps: (i) we first study the stability of the iterate DD scheme (inner solver) and estimate the condition number, and (ii) we prove the convergence of the LDD scheme (outer solver), show the well-posedness of the discrete scheme, estimate the convergence rate and subsequently determine the optimal stabilization parameter. A key point in the analysis of the methods below are inverse inequalities.

Lemma 5.1** (Inverse inequalities).**

There exist positive constants $C_{\textnormal{dTr}},C_{\textnormal{inv}}>0$ depending only on the shape regularity of the mesh such that

[TABLE]

5.1 Analysis of the DD step

To simplify, we rewrite (4.1): find $({\mathbf{u}}^{n,k}_{h,\gamma},\lambda^{n,k}_{h,\gamma})$ $\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$ so that

[TABLE]

where $\mathcal{A}_{\gamma}$ is the linearized flow system on the fracture and $s_{\gamma}$ is the flow contribution from the rock matrices

[TABLE]

The first result concerns the properties of the coupling term $s_{\gamma}$.

Lemma 5.2** (Properties of the DD operator).**

The interface bilinear form $s_{\gamma}$ is symmetric positive and semi-definite on $L^{2}(\gamma)$, and there exists a constant $C_{1}>0$ independent of $h$ such that, for all $\lambda_{h,\gamma}\in M_{h,\gamma}$,

[TABLE]

Proof.

Recalling (3.10) and taking ${\mathbf{v}}={\mathbf{u}}^{*}_{h,i}(\mu)$ and $q=p^{*}_{h,i}(\mu)$ in (3.8), $s_{\gamma}$ can be expressed as

[TABLE]

It is now easy to see that the bilinear form $s_{\gamma}$ is symmetric and positive semi-definite on $L^{2}(\gamma)$. We now show that if $s_{\gamma}(\lambda_{h,\gamma},\lambda_{h,\gamma})=0$, then $\lambda_{h,\gamma}=0$ on $M_{h,\gamma}$. Note that $s_{\gamma}(\lambda_{h,\gamma},\lambda_{h,\gamma})=0$ implies that ${\mathbf{u}}^{*}_{h,i}(\lambda_{h,\gamma})=p^{*}_{h,i}(\lambda_{h,\gamma})=0$. Again, (3.8) implies $(\mathcal{Q}_{h,i}\lambda_{h,\gamma},{\mathbf{v}}\cdot{\mathbf{n}}_{i})_{\gamma}=(\lambda_{h,\gamma},{\mathbf{v}}\cdot{\mathbf{n}}_{i})_{\gamma}=0$ for any ${\mathbf{v}}\in\mathbf{V}_{h,i}$. Thus, we can find some ${\mathbf{v}}$ so that ${\mathbf{v}}\cdot{\mathbf{n}}_{i}=\mathcal{Q}_{h,i}\lambda_{h,\gamma}$ and then $||\mathcal{Q}_{h,i}\lambda_{h,\gamma}||_{\gamma}=0$. Finally, (3.5) shows that $\lambda_{h,\gamma}=0$ on $\gamma$. We now infer the upper bound on $s_{\gamma}$. The assumption (A2) directly implies

[TABLE]

The definition (3.10) of $s_{\gamma}$ gives

[TABLE]

This result together with (5.6) leads to the upper bound in (5.5). We prove the lower bound by induction. We let $(\psi_{i},{\mathbf{r}}_{i})$, $i\in\{1,2\}$, be the solution of the auxiliary subdomain problem

[TABLE]

For fracture network with immersed fractures or for subdomains with $\Gamma_{i}=\emptyset$, $\lambda_{h,\gamma}$ approximates the pressure on $\gamma$ determined up to a constant. This constant is fixed by a zero mean value constraint for $M_{h,\gamma}$ [25, 9]. Thus, the auxiliary problem is well-posed since $({\mathbf{r}}_{i}\cdot{\mathbf{n}}_{i},1)_{\partial\Omega_{i}}=(\mathcal{Q}_{h,i}\lambda_{h,\gamma},1)_{\partial\Omega_{i}}=0$. Now, we choose ${\mathbf{v}}=\tilde{\Pi}_{i}{\mathbf{r}}_{i}$ in (3.8), to obtain

[TABLE]

where we used (5.7), assumption (A2) and the elliptic regularity (3.4) and $||{\mathbf{r}}_{i}||_{1/2,\Omega_{i}}\lesssim||\mathcal{Q}_{h,i}\lambda_{h,\gamma}||_{\gamma}$. The bound (5.9) in combination with (5.6)-(5.7) and (3.5) delivers the lower bound in (5.5).

It is interesting to study the robustness of Algorithm 4.1 and Algorithm 4.3 for the limiting case $\alpha_{\gamma}\rightarrow\infty$ in the transmission conditions (2.3), which is corresponding to a continuous pressure over the fracture interface.

Lemma 5.3** (Parameter-robustness

($\alpha_{\gamma}\rightarrow\infty$)).**

In the case of continuous pressure across $\gamma$, there exists a constant $C_{2}>0$ such that, for all $\lambda_{h,\gamma}\in M_{h,\gamma}$,

[TABLE]

Proof.

Recalling the definition (3.10) of $s_{\gamma}$ and using (5.1), we have

[TABLE]

This result together with (5.6) and (5.7) leads to the upper bound in (5.10). By inspection of the proof of Lemma 5.2, starting as in (5.9) we promptly get the lower bound of (5.10).

In the following, we denote by $||\cdot||_{s,\gamma}$ the induced semi-norm from $s_{\gamma}$ on $L^{2}(\gamma)$,

[TABLE]

We will also consider the following discrete norms:

[TABLE]

The following estimates are obtained.

Lemma 5.4** (Inverse energy estimates).**

There holds for all $({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$,

[TABLE]

Furthermore, if $\alpha_{\gamma}\rightarrow\infty$, there holds

[TABLE]

Proof.

With (5.2) and (5.5), we obtain (5.14). If $\alpha_{\gamma}\rightarrow\infty$, we make use (5.10) to get (5.15).

The following results are immediately verified.

Lemma 5.5** (Boundedness on ${A}_{\gamma}$).**

There holds for all $({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})$, $({\mathbf{v}}_{h,\gamma},\mu_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$,

[TABLE]

Lemma 5.6** (Positivity on ${A}_{\gamma}$).**

There holds for all $({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$,

[TABLE]

The above results are then used to prove the following stability estimate for $A_{\gamma}+s_{\gamma}$.

Theorem 5.7** (Stability results).**

For $({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$, we have

[TABLE]

If $\alpha_{\gamma}\rightarrow\infty$, we have

[TABLE]

Proof.

Let us first recall the inf-sup condition; given $\lambda_{h,\gamma}\in M_{h,\gamma}$, we construct ${\mathbf{r}}_{h,\gamma}\in\mathbf{V}_{h,\gamma}$ such that

[TABLE]

Let $\Psi_{\gamma}\in H^{2}_{0}(\gamma)$ satisfy $-\Delta_{\tau}\Psi_{\gamma}=\tau^{-1}\lambda_{h,\gamma}$. Take ${\mathbf{r}}_{\gamma}=-\nabla_{\tau}\Psi_{\gamma}$ and let ${\mathbf{r}}_{h,\gamma}=\Pi_{h,\gamma}{\mathbf{r}}_{\gamma}$, where $\Pi_{h,\gamma}$ is the Raviart-Thomas projection onto $\mathbf{V}_{h,\gamma}$ [36, 15]. Then $\nabla_{\tau}\cdot{\mathbf{r}}_{h,\gamma}=\Pi_{h,\gamma}\nabla_{\tau}\cdot{\mathbf{r}}_{\gamma}=\tau^{-1}\lambda_{h,\gamma}$ and $b_{\gamma}({\mathbf{r}}_{h,\gamma},\lambda_{h,\gamma})=||\tau^{-\frac{1}{2}}\lambda_{h,\gamma}||_{\gamma}^{2}$. Finally, $||{\mathbf{r}}_{h,\gamma}||_{\gamma}^{2}\leq C||{\mathbf{r}}_{\gamma}||_{1,\gamma}^{2}\leq C||\Psi_{\gamma}||_{2,\gamma}^{2}\leq C(\gamma)||\tau^{-\frac{1}{2}}\lambda_{h,\gamma}||_{\gamma}^{2}$. Set $\delta_{1},\delta_{2}>0$, and let ${\mathbf{v}}_{h,\gamma}={\mathbf{u}}_{h,\gamma}-\delta_{2}{\mathbf{r}}_{h,\gamma}$ and $\mu_{h,\gamma}=\lambda_{h,\gamma}+\delta_{1}\tau\nabla_{\tau}\cdot{\mathbf{u}}_{h,\gamma}$, where ${\mathbf{r}}_{h,\gamma}$ from (5.20). We get

[TABLE]

For the first term on the right-hand side of (5.21), we obtain using estimate (5.17) together with (5.12),

[TABLE]

For the second term, we get for all $\epsilon_{1}>0$,

[TABLE]

where we have used the continuity of $s_{\gamma}$, i.e., $s_{\gamma}(\lambda_{h,\gamma},\mu_{h,\gamma})\leq||\lambda_{h,\gamma}||_{s,\gamma}||\mu_{h,\gamma}||_{s,\gamma}\leq\alpha_{\gamma}||\lambda_{h,\gamma}||_{\gamma}||\mu_{h,\gamma}||_{\gamma}.$

For the last term, using (A2) together with (5.20) (first equation), we obtain for all $\epsilon_{2}>0$,

[TABLE]

Thus, with (5.20) (second equation),

[TABLE]

Collecting the previous results we get

[TABLE]

Now, let us choose the parameters $\epsilon_{i}$ and $\delta_{i}$ such that all the norms in the previous inequality are multiplied by positive coefficients. We choose $\epsilon_{1}=1/(1+\tau\alpha_{\gamma})$ and $\delta_{1}=1/(1+\tau\alpha_{\gamma})^{2}$, and then $\epsilon_{2}={2}/[{C(\gamma)(C_{\mathbf{K},\gamma}+L_{\gamma})}]$ and $\delta_{2}={2}/[{C(\gamma)(C_{\mathbf{K},\gamma}+L_{\gamma})+1}]$, to get

[TABLE]

Thus,

[TABLE]

From the choice of ${\mathbf{v}}_{h,\gamma}$ and $\mu_{h,\gamma}$, we have that

[TABLE]

With simple calculations, it is inferred that

[TABLE]

This implies that we have $||({\mathbf{v}}_{h,\gamma},\mu_{h,\gamma})||_{1,\tau,\star}\leq 3||({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})||_{1,\tau,\star}$. This together with (5.22) leads to (5.18). For $\alpha_{\gamma}\rightarrow\infty$, we repeat the same proof while using (5.10) instead of (5.5) to get (5.19).

Lemma 5.8** (Well-posedness of the DD scheme).**

*The domain decomposition scheme (5.3) is well-posed, and all eigenvalues of the induced system $\mathcal{A}_{\gamma}+s_{\gamma}$ are bounded away from zero. *

Proof.

This directly follows from non-singularity of $\mathcal{A}_{\gamma}+s_{\gamma}$ and estimates in Theorem 5.7.

Let us comment on the robustness of the stability estimate in Theorem 5.7. First, (5.18) states that, regardless of the choice of the space and time discretization, the stability constant with respect to $||({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})||_{1,\tau,\star}$ is independent of the coefficients $\mathbf{K}$, $\mathbf{K}_{\gamma}$, and the stabilization parameter $L_{\gamma}$. One can also show that this estimate is asymptotically robust and bounded independently of $(\tau,\alpha_{\gamma},h)\rightarrow 0$ and the stability constant tends to $1/6$. The only issue can happen having a large coefficient $\alpha_{\gamma}$, but this case is resolved in (5.19). Therein, as the ratio ${\tau}/{h}\rightarrow 0$, the stability constant is approximately $1/6$.

Following the approach of Ern and Guermond [22], we now provide an estimate for the condition number of the stiffness matrix associated with the domain decomposition system $\mathcal{A}_{\gamma}+s_{\gamma}$. This condition number estimate is important in our analysis as any algorithm is stable if every step is well-conditioned. This will also encourage the development of the flux basis framework in Section 7. Let us first introduce some basic notation in order to provide the definition of the condition number. We recall the stiffness matrix $\mathcal{A}_{\textnormal{DD}}$ introduced in (4.2) associated with the domain decomposition scheme (5.3),

[TABLE]

for all $({\mathbf{u}}_{h,\gamma},\lambda_{h,\gamma})$, $({\mathbf{v}}_{h,\gamma},\mu_{h,\gamma})\in\mathbf{V}_{h,\gamma}\times M_{h,\gamma}$, where $(V,W)_{N}\vcentcolon=\sum_{i=1}^{N}V_{i}W_{i}$ denotes the inner product in ${\amsmathbb R}^{N}$ and $|V|_{N}^{2}\vcentcolon=(V,V)_{N}$ is the corresponding Euclidean norm. The condition number is defined by $\kappa(\mathcal{A}_{\textnormal{DD}})\vcentcolon=|\mathcal{A}_{\textnormal{DD}}|_{N}|\mathcal{A}_{\textnormal{DD}}|^{-1}_{N}$ where

[TABLE]

We recall the following estimate that holds true for a conforming, quasi-uniform mesh $\mathcal{T}_{h}$ [22]; there exists $c_{\mu},C_{\mu}>0$ such that the following equivalence holds

[TABLE]

Theorem 5.9** (Condition number estimate).**

​ The condition number $\kappa(\mathcal{A}_{\textnormal{DD}})$ of (5.3) is bounded as

[TABLE]

Furthermore, if $\alpha_{\gamma}\rightarrow\infty$,

[TABLE]

Proof.

By definition (5.24), using (5.14) and (5.25), we have for all $V,W\in{\amsmathbb R}^{N}$,

[TABLE]

Consequently, $|\mathcal{A}_{\textnormal{DD}}|_{N}\lesssim\max((1+C_{\textnormal{inv}}c_{\mathbf{K},\gamma}\tau h^{-2}),(1+\alpha_{\gamma}\tau))h^{d}.$ On the other hand,

[TABLE]

hence $|V|_{N}\lesssim 6(1+\tau\alpha_{\gamma})^{2}h^{-d}|\mathcal{A}_{\textnormal{DD}}V|_{N}$. Setting $V=\mathcal{A}_{\textnormal{DD}}^{-1}W$, we conclude that $|\mathcal{A}_{\textnormal{DD}}^{-1}|_{N}\lesssim 6(1+\tau\alpha_{\gamma})^{2}h^{-d}$. Combining estimates for $|\mathcal{A}_{\textnormal{DD}}^{-1}|_{N}$ and $|\mathcal{A}_{\textnormal{DD}}|_{N}$ we get (5.26). The estimate (5.27) for $\alpha_{\gamma}\to\infty$ is obtained similarly by using (5.15) and (5.19).

5.2 Convergence of MoLDD-scheme

The second step of our analysis is to prove the convergence of Algorithm 4.1. The idea is to show that this algorithm is a contraction and then apply the Banach fixed-point theorem [44]. For that purpose, define $\delta_{{\mathbf{u}},h}^{k}={\mathbf{u}}^{n,k}_{h,\gamma}-{\mathbf{u}}^{n,k-1}_{h,\gamma}$ and $\delta_{\lambda,h}^{k}=\lambda^{n,k}_{h,\gamma}-\lambda^{n,k-1}_{h,\gamma}$ .

Theorem 5.10** (Convergence of MoLDD-scheme).**

Assuming that Assumptions (A1)–(A5) hold true and that $L_{\gamma}(\zeta)={L_{\xi}}/{2(1-\zeta)}$, with a parameter $\zeta\in[0,1)$, Algorithm 4.1 defines a contraction given by

[TABLE]

*and the limit is the unique solution of (3.11).

Proof.

By subtracting (4.1) at step $k$ from the ones at step $k-1$, we obtain

[TABLE]

Taking ${\mathbf{v}}=\tau\delta_{{\mathbf{u}},h}^{k}$ in (5.29a), $\mu=\delta_{\lambda,h}^{k}$ in (5.29b) and summing the equations, we get

[TABLE]

By choosing $L_{\gamma}={L_{\xi}}/{2(1-\zeta)}$, we immediately obtain (5.28). The inequality (5.28) implies that the sequence $\delta_{\lambda,h}^{n,k}{{\rightarrow 0}}$ in $L^{2}(\gamma)$ and $\delta_{{\mathbf{u}},h}^{k}{{\rightarrow 0}}$ in $\mathbf{L}^{2}(\gamma)$. Now we choose $\mu=\nabla_{\tau}\cdot\delta_{{\mathbf{u}},h}^{k}$ in (5.29b) to obtain

[TABLE]

Thus,

[TABLE]

Hence, by (5.28), we have $||\nabla_{\tau}\cdot\delta_{{\mathbf{u}},h}^{k}||_{\gamma}$ tends to [math] in $L^{2}(\gamma)$. This shows that $\delta_{{\mathbf{u}},h}^{k}$ tends to [math] in $\mathbf{H}(\textnormal{div}_{\tau},\gamma)$.

Since (5.28) defines a contraction, by the Banach fixed-point theorem we can conclude that the sequence generated by the algorithm converges to the unique solution of the problem (3.11).

Corollary 5.11** (Optimal MoLDD-convergence rate).**

If $\xi_{m}>0$, the minimum of the convergence rate of Algorithm 4.1 is reached for the optimal parameter

[TABLE]

Proof.

Plugging $L_{\gamma}={L_{\xi}}/{2(1-\zeta)}$ in the contraction estimate (5.28) leads to

[TABLE]

which clearly can be minimal when choosing the optimal value of $\zeta$. To calculate $\zeta^{*}$, we differentiate (5.33) with respect to $\zeta$ and infer the resulting roots and we find that the minimum of (5.33) is obtained for the optimal choice given by (5.32a). Replacing back the resulting value into $L_{\gamma}(\zeta)$ gives (5.32c).

Lemma 5.12** (Well-posedness of the mixed-dimensional problem).**

*There exists a unique solution to the mixed-dimensional problem (3.6). *

Proof.

Problem (3.11) is equivalent to (3.6). Since we know from Theorem 5.10 that (3.11) has a unique solution, this equivalence implies that (3.6) is uniquely solvable.

The rate of convergence (5.32b) depends only on the strength of the non-linearity by means of the Lipschitz constant $L_{\xi}$, the lower bound $\xi_{m}$ and the fracture permeability $\mathbf{K}_{\gamma}$. Particularly, the rate is independent of the fracture-matrix coupling parameter $\alpha_{\gamma}$, the mesh size $h$ and the time step $\tau$. Moreover,

the convergence of MoLDD-scheme is global, i.e. independent of the initialization and particularly of the inner DD solver. Nevertheless, it is beneficial if one starts MoLDD-scheme iterations with the solution of the last time step.

6 Analysis of ItLDD-scheme

We turn now to the analysis of the iterative LDD-scheme (Algorithm 4.3). In contrast to MoLDD-scheme, in which two levels of calculations (Linearization+DD) are necessary to achieve the required solution, the iterative LDD-scheme treats simultaneously the non-linearity and DD. Thus, the next result is to be understood as the convergence for the combined Linearization-DD processes. We again denote $\delta_{{\mathbf{u}},h}^{k}\vcentcolon={\mathbf{u}}^{n,k}_{h,\gamma}-{\mathbf{u}}^{n}_{h,\gamma}$ and $\delta_{\lambda,h}^{k}\vcentcolon=\lambda^{n,k}_{h,\gamma}-\lambda^{n}_{h,\gamma}$.

Theorem 6.1** (Convergence of ItLDD-scheme).**

Assuming that Assumptions (A1)–(A5) hold true and that $L_{\gamma,u}(\zeta)={L_{\xi}}/{2(1-\zeta)}$, where $\zeta$ is a parameter to be optimized in $[0,1)$, and $L_{\gamma,p}\geq\alpha_{\gamma}$, the ItLDD-scheme given by Algorithm 4.3 is linearly convergent. There holds

[TABLE]

Proof.

By subtracting (4.1) at the iteration $k$ from (3.11), we obtain

[TABLE]

Taking ${\mathbf{v}}=\tau\delta_{{\mathbf{u}},h}^{k}$ in (6.2a) and $\mu=\delta_{\lambda,h}^{k}$ in (6.2b), and summing up the equations gives

[TABLE]

For any $\zeta\in[0,1)$, this is equivalent to

[TABLE]

We apply the lower bound in the last term of the left-hand side and then use the monotonicity and Lipschitz continuity of the operator $\xi$, followed by Cauchy-Schwarz and Young’s inequalities in the second term of the right-hand side, to get

[TABLE]

the continuity of $s_{\gamma}$ we get

[TABLE]

Applying Young’s inequality to (6.6), plugging the result in (6) and choosing $L_{\gamma,u}={L_{\xi}}/{2(1-\zeta)}$ gives

[TABLE]

We let $L_{\gamma,p}\geq\alpha_{\gamma}$, to obtain the estimate (6.1) . We repeat the same techniques as in (5.31), to get that $||\nabla_{\tau}\cdot\delta_{{\mathbf{u}},h}^{k}||_{\gamma}$ tends to [math] in $L^{2}(\gamma)$. This shows that $\delta_{\lambda,h}^{n,k}$ tends to 0 in $L^{2}(\gamma)$ and $\delta_{{\mathbf{u}},h}^{k}$ tends to [math] in $\mathbf{H}(\textnormal{div}_{\tau},\gamma)$.

Our contraction estimate shows that the strength of the non-linearity and the matrix-fracture (DD) coupling controls the convergence rate. In practice, the contraction factor is better if we take into account the energy norm ${\tau}||\delta_{\lambda,h}^{k}||_{s,\gamma}^{2}/2$ using the bound (5.5). Since we assume $L_{\gamma,p}\geq\alpha_{\gamma}$, we have to study the robustness of the algorithm when $\alpha_{\gamma}\rightarrow\infty$.

Lemma 6.2** (Contraction-robustness).**

*Assuming continuous pressure across $\gamma$ ($\alpha_{\gamma}\rightarrow\infty$), then let $L_{\gamma,u}(\zeta)={L_{\xi}}/{2(1-\zeta)}$ with $\zeta$ to be chosen in $[0,1)$, and $L_{\gamma,p}\geq{C_{\textnormal{dTr}}^{2}}/({c_{\mathbf{K}}}h)$, the contraction (6.1) holds true for the ItLDD-scheme in Algorithm 4.3. *

Proof.

Recall the estimate (6) which holds true in that case. We then estimate the coupling term $|s_{\gamma}(\delta_{\lambda,h}^{k}-\delta_{\lambda,h}^{k-1},\delta_{\lambda,h}^{k})|$ with the help of (5.10),

[TABLE]

We apply Young’s inequality to (6.6) and replace the result in (6), while choosing $L_{\gamma}={L_{\xi}}/{2(1-\zeta)}$,

[TABLE]

Choosing $L_{\gamma,p}\geq{C_{\textnormal{dTr}}^{2}}/({c_{\mathbf{K}}}h)$ gives the contraction (6.1). The rest of the proof is as in Theorem 6.1.

Finally, we give alternative convergence results when $(h,{1}/{\alpha_{\gamma}})\rightarrow 0$, leading to extremely large stabilization parameter $L_{\gamma,p}$ which deteriorates the convergence rate of ItLDD scheme. These results are important to show the robustness of the ItLDD scheme for extreme physical or discretization situations.

Proposition 6.3** (Alternative convergence results).**

If $L_{\gamma,p}=0$, and $L_{\gamma,u}={L_{\xi}}/{2(1-\zeta)}$ with $\zeta\in[0,1)$, Algorithm 4.3 is convergent under the constraint on the time step $\tau\leq{1}/{\alpha_{\gamma}}$. The following estimate for Algorithm 4.3 holds true and defines a contraction

[TABLE]

Moreover, if $\alpha_{\gamma}\rightarrow\infty$, Algorithm 4.3 is convergent as ${\tau}/{h}\leq{c_{\mathbf{K}}}/{C_{\textnormal{dTr}}^{2}}(=\vcentcolon{C^{-1}_{\gamma,s}})$ holds true, and the resulting estimate is a contraction given by

[TABLE]

Proof.

We let $L_{\gamma,p}=0$ in the estimate (6.3) to get

[TABLE]

With the same techniques used to get (6), we get for $L_{\gamma,u}={L_{\xi}}/{2(1-\zeta)}$ with $\zeta\in[0,1)$,

[TABLE]

The coupling term in the right-hand side is now estimated as $s_{\gamma}(\delta_{\lambda,h}^{k-1},\delta_{\lambda,h}^{k})|\leq\alpha_{\gamma}||\delta_{\lambda,h}^{k-1}||_{\gamma}||\delta_{\lambda,h}^{k}||_{\gamma}$, where we used (5.5). Applying Young inequality and inserting the result in (6.12), we get (6.10). The second estimate (6.11), when $\alpha_{\gamma}\rightarrow\infty$, is similar to (6.10), but with using (5.10) to bound the coupling term.

The constraint on ${\tau}/{h}$ is less restrictive than the constrain on $L_{\gamma,p}$ in Lemma 6.2. Moreover, the constraint $\tau\leq{1}/{\alpha_{\gamma}}$ may have the same implication on the convergence rate as taking $L_{\gamma,p}\geq\alpha_{\gamma}$ in Theorem 6.1. In practice, the choice between the two constraints may depend on the situation. All the results show a strong correlation between the Robin parameter $\alpha_{\gamma}$, the time step $\tau$ (or ${\tau}/{h}$) and the stabilization parameter $L_{\gamma,p}$.

7 Implementation of LDD schemes using Multiscale Flux Basis

Lastly, we propose an implementation of the inter-dimensional map $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$ based on the construction of a multiscale mortar flux basis [2, 28]. We want to solve the reduced scheme in Definition 3.2 for different physical and $L$-scheme parameters, for variation of non-Darcy flow problems by changing the non-linearity $\xi$, and to compare the two LDD schemes. In Section 4, we have seen that the dominant computational cost of both Algorithm 4.1 and Algorithm 4.3 comes from the subdomain solves by evaluating the action of $\mathcal{S}_{\gamma}^{\textnormal{RtN}}$ using Algorithm 4.2 (step 2(b)). We have also seen that the condition number (5.26)-(5.27) of the linearized interface problem grows with refining the grids or increasing $\alpha_{\gamma}$. Therefore, in the case of large number of iterations, we want to have an efficient method to evaluate the action of $\mathcal{S}_{\gamma}^{\textnormal{RtN}}$.

The construction of the inter-dimensional mapping is achieved by pre-computing and storing the flux subdomain responses, called multiscale flux basis (MFB), associated with each fracture pressure degree of freedom on each subdomain [29]. Let $(\Phi^{\ell}_{h,\gamma})^{\mathcal{N}_{h,\gamma}}_{\ell=1}$ be a set of basis functions on $M_{h,\gamma}$, where $\mathcal{N}_{h,\gamma}$ is the dimension of $M_{h,\gamma}$. Then, each function $\lambda_{h,\gamma}\in M_{h,\gamma}$ can be represented as $\lambda_{h,\gamma}=\sum_{\ell=1}^{\mathcal{N}_{h,\gamma}}\lambda^{\ell}_{h,\gamma}\Phi^{\ell}_{h,\gamma}$. The MFB functions corresponding to $(\Phi^{\ell}_{h,\gamma})^{\mathcal{N}_{h,\gamma}}_{\ell=1}$ are computed as follows.

Algorithm 7.1** (Assembly of the multiscale flux basis).**

Enter the basis $(\Phi^{\ell}_{h,\gamma})^{\mathcal{N}_{h,\gamma}}_{\ell=1}$. Set $\ell\vcentcolon=0$. 2. 2.

Do**

1. (a)

Increase $\ell\vcentcolon=\ell+1$. 2. (b)

Project $\Phi^{\ell}_{h,\gamma}$ on the subdomain boundary, $\lambda_{h,i}^{\ell}=\mathcal{Q}_{h,i}(\Phi^{\ell}_{h,\gamma})$. 3. (c)

Solve problem (3.8) in each subdomain $\Omega_{i}$. 4. (d)

Project the resulting flux onto the pressure space on the fracture, $\Psi_{h,\gamma,i}^{\ell}\vcentcolon=-\mathcal{Q}^{\textnormal{T}}_{h,i}{\mathbf{u}}^{*}_{h,i}(\lambda_{h,i}^{\ell})\cdot{\mathbf{n}}_{i}.$

While $\ell\leq\mathcal{N}_{h,\gamma}$.* *

Hence, the action of $\mathcal{S}^{\textnormal{RtN}}_{\gamma}$ is given by

[TABLE]

This holds true at any time step and at any iteration of MoLDD or ItLDD. The use of MFB removes the dependence between the total number of subdomain solves and the number of their iterations.

Remark 7.2** (Fracture network).**

*For the case of a fracture network, say $\gamma=\cup_{i\neq j}\gamma_{ij}$, where $\gamma_{ij}$ is the fracture between the subdomain $\Omega_{i}$ and $\Omega_{j}$, the previous basis reconstruction is then applied independently on each fracture. *

Remark 7.3** (Inner solver).**

*Our numerical examples have a single one-dimensional fracture and, thus, we only use the direct methods to solve the interface system. For large-scale problems with many fractures or in three spatial dimensions, we emphasize the need for an iterative solver, such as GMRes. However, we see from Theorem 5.9 that the condition number of the DD system depends on $h$ and $\alpha_{\gamma}$, which in turn influences the number of iterations of the iterative solver. To retain robustness, we can use a preconditioner [18, 7] or a coarse mortar space that is compensated by taking higher order mortars [50, 9]. *

8 Numerical examples

In this section, we present several test cases to show how the schemes behave (1) for different values for numerical and physical parameters (2) with coarsening/refining mortar grids (3) on extensions to other governing equations. We subsequently study the value of $L_{\gamma,opt}$ in the MoLDD scheme and the relationship between $L_{\gamma,u}$ and $L_{\gamma,p}$ in the ItLDD scheme. The performance of schemes is measured in the overall number of iterations needed for each scheme to reach the stopping criteria. In the implementation of both schemes, we consider that the solution has converged if the relative error of the fracture solution is less than $10^{-5}$, if the value at the previous iteration step is not below machine precision. Otherwise we use the absolute error. We use a direct method (LU decomposition) to solve the linearized interface problem since the size of the system is relatively small.

To keep the presentation simple, we consider the domain and several parameters in common in all the examples in relation to the first test case in [41]. The domain $\Omega\vcentcolon=(0,2)\times(0,1)$ is intersected with a fracture defined as $\gamma\vcentcolon=\{x=1\}$. On the boundaries of the rock matrix $\{x=0\}$ and $\{x=2\}$ we impose pressure boundary condition with values [math] and $1$, respectively. We set zero flux boundary condition on the rest of $\partial\Omega$. The boundary of the fracture at the tips $\{y=1\}\cap\partial\gamma$ and $\{y=0\}\cap\partial\gamma$ inherits the pressure boundary conditions from the rock matrix. The examples are set on the time interval $I=(0,1)$ with homogeneous pressure initial condition. As for the physical parameters, we take the permeability matrix for the bulk $\mathbf{K}_{i}=\mathbf{I}$, while the source terms $f_{i}$ and $f_{\gamma}$ are equal to zero.

First, we consider the Forchheimer flow model where the non linear term is $\xi({\mathbf{u}}_{\gamma})=\beta|{{\mathbf{u}}_{\gamma}}|{\mathbf{u}}_{\gamma}$. The parameter $\beta$ is a fluid dependent non-negative scalar known as the Forchheimer coefficient, and $|\cdot|$ denotes the Euclidean vector norm $|{{\mathbf{u}}_{\gamma}}|^{2}={\mathbf{u}}_{\gamma}\cdot{\mathbf{u}}_{\gamma}$. It is straightforward to see that $\xi$ is a simply increasing function and satisfies condition (A1). For more details see [30, 36] and references therein.

8.1 Stability with respect to the user-given parameters

We first study the performances of MoLDD and ItLDD solvers by varying the time step $\tau$, the mesh size $h$, and the $L$-scheme parameters $(L_{\gamma,u},L_{\gamma,p})$. We let $\mathbf{K}_{\gamma}=1$, $\alpha_{\gamma}=10^{4}$ and $\beta=1$ and according to the theoretical results, the L-scheme parameters are given by $L_{\gamma,u}\approx 1$ and $L_{\gamma,p}=10^{3}$. Results in Table 1 report the number of iterations required by the two LDD solvers while varying the mesh size $h$ and time step size $\tau$. Each column of the tables represent results for a time step $n$.

Regardless of the choice of scheme, we can observe that the number of iterations is independent from the mesh size and slightly dependent on the time step size. The reason for the latter might be related to the fact that we consider the solution at previous iteration as the initial guess for the next iteration. Thus, by decreasing the time step size, the variation of the solution between steps varies less and so the number of iterations. Overall, the sequential ItLDD and the monolithic MoLDD solvers behave similarly; one can also see a slightly better results for the iterative solver in Table 1 (left). Note that any comparison of the two solvers does not make sense for the simple reason that the amounts of stabilization fixed by $L_{\gamma,p}$ and $L_{\gamma,u}$ are not yet optimal. Another explanation is that the amount of stabilization in the monolithic solver MoLDD is set solely by $L_{\gamma,u}$, in contrast to the iterative solver ItLDD where $L_{\gamma,p}$ and $L_{\gamma,u}$ are used.

Finally, we recall that with the use of the MFB, the computational costs of the two solvers is practically the same. The main cost is done offline, which is mostly related to the number of mortar degrees of freedom. As an example, the computational cost needed to draw the results for $h=2^{-5}$ of the two right tables in Table 1 is approximately equal to 96 subdomain solves $(\textnormal{Num. of DOF}*\textnormal{Num. of subdo.}+2N)$, where two solves per time step are required to form the right-hand side in (4.1) (for MoLDD) and (4.3) (for ItLDD). Without MFB, the cost should be $\sum_{n=1}^{N}\sum_{k=1}^{N^{n}_{\textnormal{Lin}}}N_{\textnormal{dd}}^{k}+2N$, where $N^{n}_{\textnormal{Lin}}$ is the number of iterations of the $L$-scheme, and $N_{\textnormal{dd}}^{k}$ denotes the number of DD iterations . Thus, if we assume a fixed $N_{\textnormal{dd}}^{k}$ along all the simulation, say $N_{\textnormal{dd}}^{k}=2$, this number will be at least 1012 subdomain solves, so that with MFB we make a save of approximately $91\%$ of the total subdomain solves.

In Figure 1a, we plot the number of iterations of the user-given $L_{\gamma,u}$ in MoLDD solver. We consider 100 values of $L_{\gamma,u}$, from [math] to $2.5$ with uniform step $0.025$. The other parameters are fixed as follows, $\beta=1$, $h=0.125$ and $\tau=0.2$. The graph in this figure behaves very similarly to what is usually observed for the $L$-type schemes (a typical V-shape graph), highlighting a numerically optimal value $L_{\gamma,opt}$ between $0.5$ and $1$. The number of iterations increases for small and large values of $L_{\gamma,u}$. This behavior is common for all time steps. We expect such a behavior when choosing $L_{\gamma,u}$ close to zero because it directly influences the contraction factor in (5.28). We can also see that the identified parameter $L_{\gamma,opt}\approx 1$ is close to the optimal one.

In Figure 1b, we show the performance of the ItLDD solver with regards to changing parameters $L_{\gamma,u}$ and $L_{\gamma,p}$. We consider $L_{\gamma,u}$ taking 50 values uniformly distributed on the interval $(0,2.5)$, while $L_{\gamma,p}=10^{x}$, where $x$ are 21 equidistant values on the interval $(2.2,4.2)$ with step $0.1$. We can observe on the surface plots that there is a global minimum that determines the optimal choice for $L_{\gamma,u}$ and $L_{\gamma,p}$. For example, the minimum number of iterations for this flow model is $5$ for $L_{\gamma,u}$ between $0.59$ and $1.1$ and $\log(L_{\gamma,p})$ between $2.8$ and $3$, in all time steps. Similar to MoLDD, the number of iterations required by the ItLDD solver increases when the L-scheme parameters assume low values. Particularly, the scheme diverges when $L_{\gamma,p}\leq 10^{2}$. In the analysis of the scheme we require that $L_{\gamma,p}\geq\alpha_{\gamma}$, but the lower values also allow a good convergence behavior, concluding that the theoretical lower bound is possibly too strict, but it certainly exists. Therefore, in practice, we can slightly relax the bounds on the L-scheme parameters to still obtain good performance of the solver. It is also relevant to mention that the normal permeability constant $\alpha_{\gamma}=10^{4}$ is sufficiently large to apply the limit case results in 6.2.

Crucially, the computational cost needed to draw Figure 1a and 1b, is exactly equal to only one realization with fixed $(L_{\gamma,u},L_{\gamma,p})$, confirming the utility of the MFB on fixing the total cost and avoiding any computational overhead if these parameters are not optimal.

8.2 Robustness with respect to the physical

parameters

We want to show now the robustness of the algorithms with respect to $\alpha_{\gamma}$ and $\beta$; $\alpha_{\gamma}$ controls the strength of the fracture-matrix coupling, while $\beta$ controls the strength of the non-linearity. We fix the mesh size $h=0.125$ and the time step $\tau=2^{-3}$.

In Table 2 (top), we study the dependency of the number of iterations on the Forchheimer coefficient $\beta$. The LDD solvers show a weak dependency of the number of iterations on the values of $\beta$, giving slightly better results for larger values. Overall, the MoLDD solver performs slightly better then the ItLDD. Bear in mind that changing $\beta$ directly influences $L_{\gamma,u}$. This shows that this parameter should be optimized in accordance to the given value of $\beta$. Again, we suggest that the decrease in number of iterations over time steps may be due to using the previous iteration solution as the initial guess in the subsequent iteration. Moreover, all the simulations in Table 2 (top) are run with a fixed computational cost, due to using MFB. Thus, strengthening or changing the non-linearity effects that may increase the number of iterations if $L_{\gamma,u}$ or $L_{\gamma,p}$ are not carefully set, has no practical effects on the total computational cost. We can conclude that the two solvers remain robust when strengthening the non-linearity effects.

In Table 2 (bottom), we show the dependency of the number of iterations on the Robin parameter $\alpha_{\gamma}$. Clearly, the number of iterations remains stable when strengthening or weakening matrix-fracture coupling, confirming and concluding the robustness of both schemes with respect to $\alpha_{\gamma}$. A solution is reported in Figure 2. The computational cost in Table 2 (bottom), any change of $\alpha_{\gamma}$ requires re-computing the MFB. However, this cost remains fixed when running and comparing the two LDD solvers for a fixed $\alpha_{\gamma}$.

8.3 Flexibility of coarsening/refining the mortar grids

In this set of simulations, we consider the case of weak inter-dimensional coupling by fixing $\alpha_{\gamma}=1$, with a low permeable fracture with $\mathbf{K}_{\gamma}=10^{-4}{\textbf{I}}$. We fix the following parameters: $h=2^{-5}$ (on the matrix), $L_{\gamma,u}=1$, $L_{\gamma,p}=2\cdot 10^{2}$ and $\beta=1$. We allow for a coarse scale of the mortar grids on the fracture; $h_{\gamma}=2^{-3},\,h_{\gamma}=2^{-4},\,h_{\gamma}=2^{-5}$, where the last choice corresponds to matching grids on the fracture. In Table 3, we plot the resulting number of iterations required by each LDD solver. Particularly, we can see that the sequential ItLDD solver in the matching grids has more difficulty to converge, so the effectiveness of the MFB is more pronounced. The monolithic solver MoLDD seems to be more robust with refining the mortar grids. Here, the computational cost of the construction of the inter-dimensional operator benefits from fewer mortar degrees on the fracture.

Example of a solution is depicted in Figure 3, where we can see that conforming and non-conforming grids (with $h_{\gamma}=2^{-3}$) on the fracture give indistinguishable results.

8.4 Extension to other flow models: the Cross model

The aim of this test case is to show that our LDD solvers can be applied to more general flow models. On the fracture, we assume the Cross flow model to relate $p_{\gamma}$ and ${\mathbf{u}}_{\gamma}$. We have the non-linear term given by

[TABLE]

The parameters $0\leq\omega_{\infty}<\omega_{0}$, ${{c_{\omega}}}$ and $r$ are positive scalars related to the rheology of the considered liquid. In (2.2a), $\mathbf{K}_{\gamma}$ is now replaced by $\omega_{\infty}$. We let $\omega\vcentcolon=\omega_{\infty}-\omega_{0}$ and set $\omega_{0}=2$, $\omega_{\infty}=1$, ${{c_{\omega}}}=1$, and $r=1.5$. It is possible to verify that $\xi$ satisfies the assumption (A1). See [24, 23] and the references therein.

We choose the iterative solver ItLDD and re-compute the simulations of Subsection 8.1 and 8.2 for the Cross flow model. We set then $L_{\gamma,u}=L_{\xi}/2=0.5$ and $L_{\gamma,p}=\alpha_{\gamma}=10^{3}$ as derived from the theory. The results (not shown) demonstrate first the stability of the ItLDD solver with respect to the parameters $h$ and $\tau$. Crucially, all the simulations in this example do not require additional computational cost (except fracture solves), as we use the same MFB inherited from the Forchheimer model. We set $h=2^{-5}$ with slightly coarse grids on the fracture $h_{\gamma}=2^{-4}$ and $\tau=2^{-4}$.

In Figure 4a, we show the results for the ItLDD solver on a set of realizations of $(L_{\gamma,u},L_{\gamma,p})$. The results do not differ greatly comparing to the case of Forchheimer’s flow model. The convexity of the surface plots in all time steps is clear giving away an optimal combination of values for $L_{\gamma,u}$ and $L_{\gamma,p}$. For example, we can find minimum of $5$ iterations for $L_{\gamma,u}$ between approximately $0.73$ and $1.57$, and $L_{\gamma,p}$ between $10^{2.8}$ and $10^{3}$. Note that the parameters prescribed by the theoretical results, $L_{\gamma,u}=0.5$ and $L_{\gamma,p}=10^{3}$, form a good candidate in this simulation. Finally, we mention that we can use an optimization process, as detailed in [46], in order to get the optimal values. Precisely, the fact that the choice of the stabilization parameters is independent of the mesh size, one can then run the LDD solver on a coarse spatial mesh and one time step, and study the stabilization parameters in specific intervals centred around the theoretical values. The parameters that give the lowest number of iterations are then used for the real computations. This “brute-force” optimization is simple to do in practice when using the MFB.

In Table 4, we consider to test the dependency of the number of iterations on the rheology parameters of the flow model. We provide results of several tests on $\omega$, ${{c_{\omega}}}$, and $r$. While testing for one of the parameters, the other two are fixed to either $\omega=1$, ${{c_{\omega}}}=1$ or $r=1.5$. We can observe that $\omega$ strongly influences the performance of both methods making it difficult to converge when $\omega$ gets larger, that is, when the non-linearity is stronger. For larger values of $\omega$ the number of iterations increases drastically, suggesting the necessity to adjust the $L$-scheme parameters as well as to use the MFB. The number of iterations was less dependent of the parameter ${{c_{\omega}}}$. This parameter itself contributes less to the strength of the non-linearity in comparison to $\omega$, and, thus, influencing less the performance of the solver. Finally, we can again notice a moderate dependency of number of iterations on parameter $r$. This is especially shown when $r>2$ and the exponent on the vector norm of ${\mathbf{u}}_{\gamma}$ becomes negative. Thus, the non-linear flow function $\xi$ is exponential in the values of ${\mathbf{u}}_{\gamma}$ and accounts for the very fast flow in the fractures. We finally recall that the robustness study drawn in Table 4 has the cost of one realization with fixed-parameters, confirming the role of the MFB in our solvers. For the robustness of LDD solvers with respect to the matrix-fracture coupling effects induced by the parameter $\alpha_{\gamma}$, we have seen that both solvers are robust when strengthening or weakening the coupling effects (results not shown). Example of solution is reported in Figure 4b.

9 Conclusions

In this study, we have presented two new strategies to solve a compressible single-phase flow problem in a porous medium with a fracture. In the porous medium, we have considered the classical Darcy relation between the velocity and the pressure while, in the fracture, a general non-linear law. We employ the L-scheme to handle the non-linearity term, but also to treat the inter-dimensional coupling in the second proposed algorithm. To further achieve computational speed-up, the linear Robin-to-Neumann co-dimensional map is constructed in an offline phase resulting in a problem reduced only to the fracture interface. This approach allows to change the fracture parameters, or the fracture flow model in general, without the need to recompute the problem associated with the rock matrix. We have shown the existence of optimal values for the L-scheme parameters, which are validated through several numerical tests. Future developments can be explored towards domain decomposition in time, where fast and slow fractures are solved asynchronously.

Acknowledgements

We acknowledge the financial support from the Research Council of Norway for the TheMSES (no. 250223) and ANIGMA projects (no. 244129/E20) through the ENERGIX program.