Decoupling mixed finite elements on hierarchical triangular grids for parabolic problems

TL;DR

This paper introduces a parallelizable mixed finite element method on hierarchical triangular grids for parabolic flow problems, combining domain decomposition, Raviart-Thomas elements, and Crank-Nicolson time integration.

Contribution

It presents a novel non-overlapping domain decomposition approach with cell-centered pressures and parallelizable subdomain systems for efficient parabolic problem solving.

Findings

Efficient parallel solution of flow problems.

Cell-centered pressure scheme with a 10-point stencil.

Stable and accurate time integration with Crank-Nicolson.

Abstract

In this paper, we propose a numerical method for the solution of time-dependent flow problems in mixed form. Such problems can be efficiently approximated on hierarchical grids, obtained from an unstructured coarse triangulation by using a regular refinement process inside each of the initial coarse elements. If these elements are considered as subdomains, we can formulate a non-overlapping domain decomposition method based on the lowest-order Raviart-Thomas elements, properly enhanced with Lagrange multipliers on the boundaries of each subdomain (excluding the Dirichlet edges). A suitable choice of mixed finite element spaces and quadrature rules yields a cell-centered scheme for the pressures with a local 10-point stencil. The resulting system of differential-algebraic equations is integrated in time by the Crank-Nicolson method, which is known to be a stiffly accurate scheme. As a result, we obtain independent subdomain linear systems that can be solved in parallel. The behaviour of the algorithm is illustrated on a variety of numerical experiments.

Tables

Table 1. Table 1: Pressure errors in ℓ ∞ ​ ( ℓ 2 ) superscript ℓ superscript ℓ 2 \ell^{\infty}(\ell^{2}) -norm obtained for a nested collection of spatial meshes 𝒯 h ℓ superscript subscript 𝒯 ℎ ℓ \mathcal{T}_{h}^{\ell} , with ℓ = 3 , 4 , 5 , 6 ℓ 3 4 5 6 \ell=3,4,5,6 , and various decreasing time steps τ 𝜏 \tau .

${\Vert r_{h}p-P_h\Vert }_{{\ell }^{\infty }\left({\ell }^2\right)}$	$\ell =3$	$\ell =4$	$\ell =5$	$\ell =6$
$\tau =4.0$E-1	1.3948E-2	8.7120E-3	7.5778E-3	7.7539E-3
$\tau =2.0$E-1	9.9587E-3	2.9405E-3	2.1206E-3	1.9511E-3
$\tau =1.0$E-1	1.0282E-2	2.7437E-3	7.5382E-4	4.2576E-4
$\tau =5.0$E-2	1.0264E-2	2.7232E-3	7.0609E-4	1.8993E-4
$\tau =2.5$E-2	1.0267E-2	2.7211E-3	7.0099E-4	1.7894E-4

Table 2. Table 2: Pressure errors in ℓ ∞ ​ ( ℓ ∞ ) superscript ℓ superscript ℓ \ell^{\infty}(\ell^{\infty}) -norm obtained for a nested collection of spatial meshes 𝒯 h ℓ superscript subscript 𝒯 ℎ ℓ \mathcal{T}_{h}^{\ell} , with ℓ = 3 , 4 , 5 , 6 ℓ 3 4 5 6 \ell=3,4,5,6 , and various decreasing time steps τ 𝜏 \tau .

${\Vert r_{h}p-P_h\Vert }_{{\ell }^{\infty }\left({\ell }^{\infty }\right)}$	$\ell =3$	$\ell =4$	$\ell =5$	$\ell =6$
$\tau =4.0$E-1	3.2717E-2	2.0340E-2	1.8256E-3	1.8446E-2
$\tau =2.0$E-1	3.0892E-2	9.3021E-3	5.0352E-3	4.6947E-3
$\tau =1.0$E-1	3.2036E-2	9.4528E-3	2.6495E-3	1.0058E-3
$\tau =5.0$E-2	3.2013E-2	9.4278E-3	2.6247E-3	7.1625E-4
$\tau =2.5$E-2	3.2083E-2	9.4370E-3	2.6212E-3	7.1009E-4

Table 3. Table 3: Errors for the post-processed fluxes for a nested collection of spatial meshes 𝒯 h ℓ superscript subscript 𝒯 ℎ ℓ \mathcal{T}_{h}^{\ell} , with ℓ = 3 , 4 , 5 , 6 ℓ 3 4 5 6 \ell=3,4,5,6 , and various decreasing time steps τ 𝜏 \tau .

${\Vert 𝐮-𝐑\left({𝐮}_h\right)\Vert }_{{\ell }^{\infty }\left(L^2\right)}$	$\ell =3$	$\ell =4$	$\ell =5$	$\ell =6$
$\tau =4.0$E-1	5.6553E-1	2.4799E-1	1.6123E-1	1.5100E-1
$\tau =2.0$E-1	4.7129E-1	1.4528E-1	5.5718E-2	4.0916E-2
$\tau =1.0$E-1	4.9077E-1	1.4819E-1	4.4715E-2	1.4333E-2
$\tau =5.0$E-2	4.8997E-1	1.4771E-1	4.4354E-2	1.3802E-2
$\tau =2.5$E-2	4.8991E-1	1.4766E-1	4.4307E-2	1.3763E-2

Table 4. Table 4: Errors for the normal fluxes obtained for a nested collection of spatial meshes 𝒯 h ℓ superscript subscript 𝒯 ℎ ℓ \mathcal{T}_{h}^{\ell} , with ℓ = 3 , 4 , 5 , 6 ℓ 3 4 5 6 \ell=3,4,5,6 , and various decreasing time steps τ 𝜏 \tau .

${\Vert \left(𝐮-{𝐮}_h\right)\cdot 𝐧\Vert }_{{\ell }^{\infty }\left({\ell }^2\right)}$	$\ell =3$	$\ell =4$	$\ell =5$	$\ell =6$
$\tau =4.0$E-1	1.0447E-1	5.8187E-2	5.2154E-2	5.2289E-2
$\tau =2.0$E-1	8.5800E-2	2.9946E-2	1.4575E-2	1.2743E-2
$\tau =1.0$E-1	8.9694E-2	3.0881E-2	1.0513E-2	3.6249E-3
$\tau =5.0$E-2	8.9650E-2	3.0849E-2	1.0482E-2	3.5917E-3
$\tau =2.5$E-2	8.9639E-2	3.0841E-2	1.0476E-2	3.5872E-3

Table 5. Table 5: Pressure and velocity errors for the Mackinnon and Carey test example. A nested collection of spatial meshes 𝒯 h ℓ superscript subscript 𝒯 ℎ ℓ \mathcal{T}_{h}^{\ell} , with ℓ = 1 , 2 , 3 , 4 , 5 ℓ 1 2 3 4 5 \ell=1,2,3,4,5 , and a time step τ = 1 𝜏 1 \tau=1 E-01 are considered. The last row in the table contains the average orders of convergence.

	${\Vert r_{h}p-P_h\Vert }_{{\ell }^{\infty }\left({\ell }^2\right)}$	${\Vert r_{h}p-P_h\Vert }_{{\ell }^{\infty }\left({\ell }^{\infty }\right)}$	${\Vert 𝐮-𝐑\left({𝐮}_h\right)\Vert }_{{\ell }^{\infty }\left(L^2\right)}$	${\Vert \left(𝐮-{𝐮}_h\right)\cdot 𝐧\Vert }_{{\ell }^{\infty }\left({\ell }^2\right)}$
$\ell =1$	2.4218E-2	7.2240E-2	3.7393E-1	2.7329E-1
$\ell =2$	5.9113E-3	1.8738E-2	1.0901E-1	8.0644E-2
$\ell =3$	1.4637E-3	4.6708E-3	3.1850E-2	2.3434E-2
$\ell =4$	3.6745E-4	1.1656E-3	9.5547E-3	6.7958E-3
$\ell =5$	9.2331E-5	2.9124E-4	2.9657E-3	2.0113E-3
order	2.009	1.989	1.745	1.772

Equations

$$p_t+\nabla \cdot \mathbf{u}=f$$

$$\mathbf{u}=-K\nabla p$$

$$p=p_0$$

$$p=g_D$$

$$\mathbf{u}\cdot \mathbf{n}=g_N$$

$$H(\text{mbox}div;R)=\{\mathbf{v}\in (L^2(R){)}^2:\nabla \cdot \mathbf{v}\in L^2(R)\}.$$

$$\tilde{\mathbf{u}}=-G^{-1}\nabla p,$$

$$G\tilde{\mathbf{u}}=-\nabla p$$

$$\mathbf{u}=KG\tilde{\mathbf{u}}$$

$$\overline{\Omega }={\cup }_{i=1}^s{\overline{\Omega }}_i,{\Omega }_i\cap {\Omega }_j=\varnothing ,\mathrm{for}i\ne j.$$

$$V=\underset{i=1}{\overset{s}{⨁}}{V}_i,\tilde{V}=\underset{i=1}{\overset{s}{⨁}}{\tilde{V}}_i=(L^2(\Omega ){)}^2,W=\underset{i=1}{\overset{s}{⨁}}{W}_i=L^2(\Omega ).$$

$$(p_t,w)+\sum _{i=1}^s(\nabla \cdot \mathbf{u},w{)}_{{\Omega }_i}=(f,w)$$

$$(G\tilde{\mathbf{u}},\mathbf{v})=\sum _{i=1}^s\left((p,\nabla \cdot \mathbf{v}{)}_{{\Omega }_i}-⟨\lambda ,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_i}\right)-⟨g_D,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_D}$$

$$(G\mathbf{u},\tilde{\mathbf{v}})=(GKG\tilde{\mathbf{u}},\tilde{\mathbf{v}})$$

$$\sum _{i=1}^s⟨\mathbf{u}\cdot \mathbf{n},\mu {⟩}_{{\Gamma }_i}=⟨g_N,\mu {⟩}_{{\Gamma }_N}$$

$$p(0)=p_0.$$

$${\hat{V}}_h(\hat{T})=({\mathbb{P}}_0(\hat{T}){)}^2\oplus \hat{\mathbf{x}}{\mathbb{P}}_0(\hat{T}),{\hat{W}}_h(\hat{T})={\mathbb{P}}_0(\hat{T}),$$

$$w\leftrightarrow \hat{w}:$$

$$\mu \leftrightarrow \hat{\mu }:$$

$$\mathbf{v}\leftrightarrow \hat{\mathbf{v}}:\mathbf{v}=\left(\frac{1}{J_T}{B}_T\hat{\mathbf{v}}\right)\circ F_T^{-1},$$

$$V_{h,i}=\left\{\mathbf{v}\in V_i:\mathbf{v}{|}_T\leftrightarrow \hat{\mathbf{v}},\hat{\mathbf{v}}\in {\hat{V}}_h(\hat{T})\forall T\in {\mathcal{T}}_{h,i}\right\},$$

$${\tilde{V}}_{h,i}=\left\{\tilde{\mathbf{v}}\in {\tilde{V}}_i:\tilde{\mathbf{v}}{|}_T\leftrightarrow \widehat{\tilde{\mathbf{v}}},\widehat{\tilde{\mathbf{v}}}\in {\hat{V}}_h(\hat{T})\forall T\in {\mathcal{T}}_{h,i}\right\},$$

$$W_{h,i}=\left\{w\in W_i:w{|}_T\leftrightarrow \hat{w},\hat{w}\in {\hat{W}}_h(\hat{T})\forall T\in {\mathcal{T}}_{h,i}\right\}.$$

$$V_h=\underset{i=1}{\overset{s}{⨁}}{V}_{h,i},{\tilde{V}}_h=\underset{i=1}{\overset{s}{⨁}}{\tilde{V}}_{h,i},W_h=\underset{i=1}{\overset{s}{⨁}}{W}_{h,i}.$$

$$L_h=\left\{\mu \in L:\mu {|}_e\leftrightarrow \hat{\mu },\hat{\mu }\in {\hat{L}}_h(\hat{e})\forall e\subset (\partial T\cap \Gamma ),T\in {\mathcal{T}}_h\right\},$$

$$N_V=\dim (V_h)=N_e+N_e^{\Gamma }-N_e^N,$$

$$N_W=\dim (W_h)=N_T,$$

$$(p_{h,t},w)+\sum _{i=1}^s(\nabla \cdot {\mathbf{u}}_h,w{)}_{{\Omega }_i}=(f,w)$$

$$(G{\tilde{\mathbf{u}}}_h,\mathbf{v})=\sum _{i=1}^s((p_h,\nabla \cdot \mathbf{v}{)}_{{\Omega }_i}-⟨{\lambda }_h,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_i})-⟨g_D,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_D}$$

$$(G{\mathbf{u}}_h,\tilde{\mathbf{v}})=(GKG{\tilde{\mathbf{u}}}_h,\tilde{\mathbf{v}})$$

$$\sum _{i=1}^s⟨{\mathbf{u}}_h\cdot \mathbf{n},\mu {⟩}_{{\Gamma }_i}=⟨g_N,\mu {⟩}_{{\Gamma }_N}$$

$$p_h(0)={\mathcal{S}}_h{p}_0.$$

$$(G\mathbf{q},\mathbf{v}{)}_T={\left(\frac{1}{J_T}{B}_T^{T}G{B}_T\hat{\mathbf{q}},\hat{\mathbf{v}}\right)}_{\hat{T}}.$$

$$(G\mathbf{q},\mathbf{v}{)}_{Q,T}\equiv (\hat{\mathbf{q}},\hat{\mathbf{v}}{)}_{\hat{Q},\hat{T}}=\frac{|\hat{T}|}{6}\left(\sum _{i=1}^3\hat{\mathbf{q}}({\hat{\mathbf{r}}}_i)\cdot \hat{\mathbf{v}}({\hat{\mathbf{r}}}_i)+3\hat{\mathbf{q}}({\hat{\mathbf{x}}}_c)\cdot \hat{\mathbf{v}}({\hat{\mathbf{x}}}_c)\right),$$

$$({\hat{\mathbf{v}}}_i,{\hat{\mathbf{v}}}_j{)}_{\hat{Q},\hat{T}}=\{\begin{array}{ll}0,& \mathrm{if}i\ne j,\\ \frac{|\hat{T}|}{6}|e_i{|}^2,& \mathrm{if}i=j,\end{array}$$

$$(G\mathbf{q},\mathbf{v}{)}_Q=\sum _{T\in {\mathcal{T}}_h}(G\mathbf{q},\mathbf{v}{)}_{Q,T}.$$

$$(G{\tilde{\mathbf{u}}}_h,\mathbf{v}{)}_Q=\sum _{i=1}^s((p_h,\nabla \cdot \mathbf{v}{)}_{{\Omega }_i}-⟨{\lambda }_h,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_i})-⟨g_D,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_D}$$

$$(G{\mathbf{u}}_h,\mathbf{v}{)}_Q=(GKG{\tilde{\mathbf{u}}}_h,\mathbf{v})$$

$${\mathbf{u}}_h(\mathbf{x},t)=\sum _{i=1}^{N_V}{U}_{h,i}(t){\mathbf{v}}_i(\mathbf{x}),$$

Full text

1 Introduction

Let us consider the parabolic initial-boundary value problem

[TABLE]

where $\Omega\subset\mathbb{R}^{2}$ is a convex polygonal domain with boundary $\partial\Omega=\overline{\Gamma}_{D}\cup\overline{\Gamma}_{N}$ such that $\Gamma_{D}\cap\Gamma_{N}=\emptyset$. In this formulation, $K\equiv K(\mathbf{x})\in\mathbb{R}^{2\times 2}$ is a symmetric and positive definite tensor, and $\mathbf{n}$ is the outward unit vector normal to $\partial\Omega$. In the framework of flow in porous media, $p$ represents the fluid pressure, $\mathbf{u}$ denotes the Darcy velocity and $K$ is the rock permeability divided by the fluid viscosity. In this context, it is common for $K$ to be a piecewise continuous tensor, whose discontinuities are typically related to the presence of layers of different material within the flow domain.

The aim of this paper is to design a numerical scheme for approximating the solution of (1). The method will be specially conceived to handle complex geometries and discontinuous tensor coefficients while preserving the computational efficiency and ease of implementation. The use of so-called hierarchical grids permits to reach the compromise between both aspects, since they combine the flexibility of unstructured triangulations –to handle irregular boundaries, or interfaces generated by discontinuous coefficients–, with the efficiency provided by regular meshes. In order to construct a hierarchical grid, we consider an unstructured coarse partition of $\Omega$ into a certain number of triangular elements, and subsequently apply a regular refinement process to each of these coarse elements111The regular triangulation obtained inside each coarse triangle is known as a three-line mesh, since it is composed of parallel lines to the sides of such a triangle.. The use of this kind of grids has experienced a great development in recent years, due to their potential for solving a considerable class of problems in large scale simulation (see, e.g., [6, 14]). Hierarchical grids are also referred to as semi-structured grids, structured multiblock grids or patch-wise structured grids (cf. [22] and references therein).

In the sequel, the coarse elements of our hierarchy of grids will be considered as subdomains of a suitable non-overlapping decomposition of $\Omega$. For the spatial discretization, we consider an expanded mixed finite element (MFE) approximation to (1) using the lowest-order Raviart–Thomas spaces. Since these spaces are separately defined for each subdomain, we need to introduce Lagrange multiplier pressures on the subdomain interfaces (and also on the Neumann boundaries) to guarantee the continuity of the normal components of velocity vectors across subdomains. This idea, first proposed in [1, 2], will later permit to decouple the global problem into a set of independent subdomain problems that can be efficiently solved in parallel.

The key idea behind expanded methods is to introduce an auxiliary variable, called the adjusted gradient, as a third explicit unknown (see [1, 9, 13]). In this way, unlike the standard mixed method, the resulting scheme avoids inverting $K$, thus allowing for a non-negative tensor coefficient. Further, it has the additional advantage that numerical quadrature can be used to obtain a cell-centered finite difference scheme for the pressures (cf. [1, 2]). Similar ideas leading to cell-centered methods on triangular meshes can be found in [8, 11, 20]. In our case, this strategy yields a stencil-based formulation of the original problem that takes advantage of the structure of the hierarchical grid. In particular, if the material properties (given by $K$) are constant within each subdomain, a unique 10-point stencil will be associated to every interior cell in the subdomain. Moreover, this stencil will be constant regardless of the level of refinement of the hierarchical grid, thus drastically reducing the computational complexity of the problem at hand.

The expanded MFE method, combined with numerical quadrature and enhanced with Lagrange multipliers at the subdomain interfaces, achieves second-order superconvergence for the pressure variable at the centroids of the triangles (see [2]). On equilateral three-line meshes and if $K$ is the identity matrix, second-order superconvergence is also observed for the normal velocity at the midpoints of the edges (cf. [10]). As we will show experimentally, the order of convergence is slightly less than 2 for the normal velocity if a general three-line mesh and/or a general tensor $K$ are considered. The application of a post-processing technique introduced in [10] will further permit to increase the optimal first-order convergence of the velocity vector to superconvergence of order almost 2.

The continuous-in-time semidiscrete problem leads to a system of semi-explicit differential-algebraic equations (DAEs) of index 1. Specifically, the velocity and gradient can be eliminated from the semidiscrete problem via a static condensation procedure, and the resulting scheme is formulated in terms of pressure and Lagrange multiplier unknowns. This formulation comprises a system of ordinary differential equations (ODEs) subject to certain constraints, that is to say, a system of DAEs in which the pressure represents the differential component and the Lagrange multiplier is the algebraic component. In order to guarantee the same order of convergence for both components, a stiffly accurate Runge–Kutta method is used for the time integration (see, e.g., [17]). In particular, we will consider the simplest member in the Lobatto IIIA family, the well-known Crank–Nicolson method, that achieves second-order convergence for both the differential and algebraic variables (cf. [16]). At each time step, the fully discrete scheme consists of a set of linear systems that decouple across subdomains and, therefore, can be efficiently solved in parallel. This parallel solution strategy has two distinct properties: on the one hand, no iterations are required for convergence; on the other, since the use of hierarchical grids implies the same size (in terms of degrees of freedom) for each subdomain problem, a perfectly balanced workload is attained among parallel processors.

The present work extends the ideas introduced in [1, 2] for the elliptic problem to the parabolic case. Further, it provides a full derivation of the explicit expressions for the stencil coefficients, which was lacking so far. In doing so, we apply a novel strategy that relies on the three-line structure of each subdomain partition, and permits an efficient implementation of the proposed algorithm. As discussed in [22] (where such an strategy is used for Galerkin finite element discretizations of elliptic problems), the idea is based on a suitable Cartesian distribution of the subdomain nodes.

The solution of evolutionary problems of the form (1) has also been tackled in the earlier work [4]. In that paper, an overlapping domain decomposition approach for solving (1), with $g_{N}=0$, is designed and analyzed. The method combines an expanded MFE scheme on triangles with a first-order time-splitting technique222Unlike the method presented here, the scheme in [4] makes use of a domain decomposition technique that is not related to a coarse triangulation of $\Omega$. Instead, it is constructed to generate a family of partition-of-unity functions that permits to split the discrete elliptic operator into terms of lower dimension yielding a time-splitting scheme.. However, its use is restricted to three-line meshes and it is not suitable for dealing with unstructured triangulations and/or piecewise continuous tensor coefficients. The newly proposed method succeeds to circumvent these situations by enhancing the method with Lagrange multipliers; in addition, it provides second-order convergence in time.

To conclude, it is remarkable to note that this enhanced variant of the expanded MFE scheme is closely related to the so-called hybrid MFE method. The hybrid method, first proposed in [3], is defined by adding Lagrange multipliers on every single edge of the triangulation of $\Omega$ (see also [7]). By construction, the present scheme, which restricts the use of Lagrange multipliers to the subdomain interfaces, will be in general much more efficient than the hybrid technique: as long as $K$ is continuous over relatively large subdomains and/or the coarse elements in the hierarchical grid are sufficiently refined, the number of Lagrange multiplier unknowns will be much larger for the hybrid formulation (cf. [2]).

The rest of the paper is organized as follows. In Section 2, we introduce the expanded variational formulation of (1). The construction of the expanded MFE spatial discretization on hierarchical grids is discussed in Section 3. This section also includes the matrix formulation of the semidiscrete scheme leading to a cell-centered finite difference method. Section 4 is devoted to the time integration of the resulting DAE system. The numerical behaviour of the fully discrete scheme is tested on a variety of numerical experiments in Section 5. The paper ends with some concluding remarks and perspectives. In addition, Appendix A includes a full derivation of the pressure stencil associated to both an upward and a downward triangle in a three-line mesh.

2 The expanded variational form

For a domain $R\subset\mathbb{R}^{2}$, let $W^{k,p}(R)$ denote the standard Sobolev space, with $k\in\mathbb{R}$ and $1\leq p\leq\infty$. Let $H^{k}(R)$ be the Hilbert space $W^{k,2}(R)$. In the sequel, we will mainly use the spaces of square-integrable scalar and vector functions, $L^{2}(R)$ and $(L^{2}(R))^{2}$, respectively. Such spaces are endowed with the inner product $(\cdot,\cdot)_{R}$ and norm $\|\cdot\|_{R}$, so that $\|\varphi\|_{R}=(\varphi,\varphi)^{1/2}_{R}$. The subscript $R$ will be omitted whenever $R\equiv\Omega$. For a section $S$ of the domain boundary, $\langle\cdot,\cdot\rangle_{S}$ represents the $L^{2}(S)$-inner product or duality pairing. We will also use the space

[TABLE]

In order to define an expanded formulation, we need to introduce an auxiliary unknown, the so-called adjusted gradient $\tilde{\mathbf{u}}\equiv\tilde{\mathbf{u}}(\mathbf{x},t)$, given by

[TABLE]

where $G\equiv G(\mathbf{x})\in\mathbb{R}^{2\times 2}$ is a symmetric and positive definite matrix related to the geometry of $\Omega$. Considering this new variable, the equation (1b) can be rewritten as

[TABLE]

The system given by the equations (1a), (2a) and (2b), together with the corresponding initial and boundary conditions, is usually referred to as expanded mixed formulation.

Let us consider that $\Omega$ can be decomposed into a set of $s$ non-overlapping triangular subdomains $\{\Omega_{1},\Omega_{2},\ldots,\Omega_{s}\}$, defined to be shape-regular of diameter $H$ and pairwise disjoint, i.e.,

[TABLE]

Let $\partial\Omega_{i}$ denote the boundary of subdomain $\Omega_{i}$. This decomposition is assumed to be geometrically conforming in the sense that any edge of $\partial\Omega_{i}$ is either part of $\partial\Omega$ or coincides with an edge of an adjacent subdomain. If $\Omega_{i}$ and $\Omega_{j}$ are adjacent to each other, we denote $\Gamma_{ij}=\partial\Omega_{i}\cap\partial\Omega_{j}$, for $i\neq j$, assuming that $\Gamma_{ij}=\Gamma_{ji}$. Furthermore, we set $\Gamma_{i}=\partial\Omega_{i}\setminus\Gamma_{D}$ and define $\Gamma=\cup_{i=1}^{s}\Gamma_{i}$. Typically, as we will see below, these subdomains correspond to the elements in a coarse partition $\mathcal{T}_{H}$ of $\Omega$ with mesh size $H$.

Next, we define the subdomain spaces $V_{i}=H(\mathrm{div};\Omega_{i})$, $\tilde{V}_{i}=(L^{2}(\Omega_{i}))^{2}$ and $W_{i}=L^{2}(\Omega_{i})$, and the global spaces

[TABLE]

Note that, in the definition of $V$, we are relaxing the continuity of the space $H(\mathrm{div};\Omega)$ across the subdomain interfaces. Thus, we need to impose it weakly by introducing the additional space $L=H^{1/2}(\Gamma)$. The following weak form can be obtained by integrating the original equations over each $\Omega_{i}$ and summing up: Find $(\mathbf{{u}},\tilde{\mathbf{u}},p,\lambda):[0,T]\rightarrow V\times\tilde{V}\times W\times L$ such that

[TABLE]

Note that the flux continuity equation (4d) implies that $\mathbf{u}$ is, in fact, an element of $H(\mathrm{div};\Omega)$. This formulation is referred to as the macro-hybrid expanded variational form of (1) with respect to the decomposition (3) (cf. [7]).

3 The expanded mixed finite element method

In this section, we describe the spatial discretization of the variational formulation (4) on a hierarchical grid. To this end, we first define the lowest-order Raviart–Thomas spaces for each single subdomain, and then introduce a space of Lagrange multipliers that ensures the flux continuity across subdomains. The subsequent application of a suitable quadrature rule for computing certain vector inner products permits to reduce the expanded MFE method to a cell-centered finite difference scheme for the pressure, enhanced with some additional unknowns on the subdomain interfaces. A static condensation procedure yields the Schur complement form of the semidiscrete system.

3.1 A hierarchical triangular grid

Let $\mathcal{T}_{H}$ be an unstructured coarse triangulation of $\Omega$, whose elements are assumed to be the subdomains $\Omega_{i}$, for $i=1,2,\ldots,s$. This partition is constructed to adequately represent the geometry of the domain, and may also take into account physical features of the problem, such as material properties. If we divide each subdomain $\Omega_{i}$ into four congruent triangles by connecting the midpoints of its edges, a new regular mesh $\mathcal{T}_{h,i}^{1}$ is obtained per subdomain. We assume that $\mathcal{T}_{h,i}^{1}$ and $\mathcal{T}_{h,j}^{1}$ match on $\Gamma_{ij}$, for $i\neq j$, so that $\mathcal{T}_{h}^{1}=\cup_{i=1}^{s}\mathcal{T}_{h,i}^{1}$ is a conforming triangulation of $\Omega$. This regular refinement process can be subsequently applied in order to obtain a nested hierarchy of conforming meshes $\mathcal{T}_{H}\subset\mathcal{T}_{h}^{1}\subset\mathcal{T}_{h}^{2}\subset\ldots\subset\mathcal{T}_{h}^{k}\equiv\mathcal{T}_{h}$, where $\mathcal{T}_{h}=\cup_{i=1}^{s}\mathcal{T}_{h,i}$ is a mesh with the desired fine scale $h$ (cf. [22]). Here, $H=\mathrm{max}_{T\in\mathcal{T}_{H}}\mathrm{diam}(T)$ and $h=\mathrm{max}_{T\in\mathcal{T}_{h}}\mathrm{diam}(T)$, where $T$ denotes an element of either $\mathcal{T}_{H}$ or $\mathcal{T}_{h}$. Figure 1 shows an example of the coarse and fine triangulations, $\mathcal{T}_{H}$ and $\mathcal{T}_{h}\equiv\mathcal{T}_{h}^{2}$, for a polygonal domain $\Omega$. Note that, for $i=1,2,\ldots,s$, each $\mathcal{T}_{h,i}$ is a three-line mesh.

3.2 Mixed finite element spaces

Let $\hat{T}$ be the reference equilateral triangle with vertices $\hat{\mathbf{r}}_{1}=(-1,0)^{T}$, $\hat{\mathbf{r}}_{2}=(1,0)^{T}$ and $\hat{\mathbf{r}}_{3}=(0,\sqrt{3})^{T}$, and introduce a family of bijective affine mappings $\{F_{T}\}_{T\in\mathcal{T}_{h}}$ such that $F_{T}(\hat{T})=T$. We further define, for each mapping $F_{T}$, the Jacobian matrix $B_{T}$ and its determinant $J_{T}=|\det(B_{T})|$. The corresponding vertices of $T$ are denoted by $\mathbf{r}_{i}=(x_{i},y_{i})^{T}$, while the outward unit vectors normal to the edges of $\hat{T}$ and $T$ are represented by $\hat{\mathbf{n}}_{i}$ and $\mathbf{n}_{i}$, respectively, for $i=1,2,3$ (see Figure 2).

Let $\hat{V}_{h}(\hat{T})$ and $\hat{W}_{h}(\hat{T})$ be the $RT_{0}$ finite element spaces on the reference element $\hat{T}$, i.e.,

[TABLE]

where $\mathbb{P}_{0}(\hat{T})$ denotes the set of constant functions defined on $\hat{T}$. On every $\hat{e}\subset\partial\hat{T}$, we further define the space $\hat{L}_{h}(\hat{e})=\mathbb{P}_{0}(\hat{e})$. If $\hat{\mathbf{v}}\in\hat{V}_{h}(\hat{T})$ and $\hat{w}\in\hat{W}_{h}(\hat{T})$, the degrees of freedom for $\hat{V}_{h}(\hat{T})$ are chosen to be the values of $\hat{\mathbf{v}}\cdot\hat{\mathbf{n}}_{i}$ at the midpoints of the edges, for $i=1,2,3$, while that for $\hat{W}_{h}(\hat{T})$ is the value of $\hat{w}$ at the element center. Finally, the corresponding degree of freedom for $\hat{\mu}\in\hat{L}_{h}(\hat{e})$ is the value of $\hat{\mu}$ at the midpoint of the edge $\hat{e}$.

In order to transform any scalar function $\hat{w}$ on $\hat{T}$ or $\hat{\mu}$ on $\hat{e}\subset\partial\hat{T}$ to a generic element $T$ or edge $e\subset\partial T$ belonging to $\mathcal{T}_{h}$, we introduce the standard isomorphisms

[TABLE]

For vector functions $\hat{\mathbf{v}}$ on $\hat{T}$, we use the Piola transformation (cf. [23])

[TABLE]

which is defined to preserve the continuity of the normal components of velocity vectors across interelement edges. This is a necessary condition that must be fulfilled when building approximations to $H(\mathrm{div};\Omega_{i})$.

The subdomain spaces on $\mathcal{T}_{h,i}$, denoted as $V_{h,i}\times\tilde{V}_{h,i}\times W_{h,i}\subset V_{i}\times\tilde{V}_{i}\times W_{i}$, are defined to be

[TABLE]

The global spaces on $\mathcal{T}_{h}$ are thus given by

[TABLE]

Note that the vector functions in $V_{h}$ have continuous normal components on the edges between elements inside each subdomain, but the space has no continuity constraint on the edges lying on $\Gamma$. In order to recover the continuity for such edges, we introduce the space $L_{h}\subset L$ of pressure Lagrange multipliers, defined as

[TABLE]

whose elements provide an approximation to $p|_{\Gamma}$.

To determine the dimensions of the preceding spaces, we need to introduce some notations first. Let $N_{e}$ and $N_{T}$ be the number of edges and elements in $\mathcal{T}_{h}$, respectively. Moreover, let $N_{e}^{\Gamma}$, $N_{e}^{D}$ and $N_{e}^{N}$ denote the number of edges in $\mathcal{T}_{h}$ belonging to $\Gamma$, $\Gamma_{D}$ and $\Gamma_{N}$, respectively. Then, we have that

[TABLE]

Observe that, in general, $N_{e}^{\Gamma}\ll N_{e}$, that is, the number of Lagrange multipliers is much smaller than the total number of edges in $\mathcal{T}_{h}$. As we will see below, this fact makes the proposed scheme much more efficient than the class of so-called hybrid methods. For such methods, a Lagrange multiplier is introduced on every edge in $\mathcal{T}_{h}$, so that the dimension of $L_{h}$ (i.e., the number of Lagrange multiplier unknowns) is $N_{e}$ instead of $N_{e}^{\Gamma}$.

The degrees of freedom for the pressure spaces $W_{h}$ and $L_{h}$ are represented in Figure 3 for a polygonal domain $\Omega$. The unknowns $p_{h}\in W_{h}$ are located at the centers of the elements $T\in\mathcal{T}_{h}$, while the Lagrange multipliers $\lambda_{h}\in L_{h}$ live on the midpoints of the edges $e\subset(\partial T\cap\Gamma)$. These latter include the edges on the internal boundaries of the subdomains $\Omega_{i}$, for $i=1,2,\ldots,s$, and the edges on the Neumann boundaries $\Gamma_{N}$.

3.3 The semidiscrete scheme

The expanded MFE approximation to the variational form (4) reads as follows: Find $(\mathbf{{u}}_{h},\tilde{\mathbf{u}}_{h},p_{h},\lambda_{h}):[0,T]\rightarrow V_{h}\times\tilde{V}_{h}\times W_{h}\times L_{h}$ such that

[TABLE]

In this case, the flux matching condition (5d) is the discrete analogue to (4d) and guarantees that $\mathbf{u}_{h}\in H(\mathrm{div};\Omega)$. The operator $\mathcal{S}_{h}$ denotes the MFE elliptic projection as defined, e.g., in [24, 25].

Considering the left-hand side of (5c) in the previous formulation, we can define a rectangular matrix $S$ whose elements are given by $(S)_{ij}=(G\mathbf{v}_{j},\tilde{\mathbf{v}}_{i})$, where $\{\mathbf{v}_{i}\}_{i=1}^{N_{V}}$ and $\{\tilde{\mathbf{v}}_{i}\}_{i=1}^{N_{\tilde{V}}}$ are bases for $V_{h}$ and $\tilde{V}_{h}$, respectively. In turn, the left-hand side of (5b) involves the matrix $S^{T}$. Following [2], we subsequently consider $\tilde{V}_{h}=V_{h}$ with the aim of making $S$ be a square, symmetric and invertible matrix. In addition, the resulting matrix can be further diagonalized by using a suitable quadrature rule for computing the inner product $(G\mathbf{v}_{j},\tilde{\mathbf{v}}_{i})$. In this way, the method becomes a cell-centered finite difference scheme for the pressure, with as many unknowns as the number of elements in $\mathcal{T}_{h}$, $N_{T}$, enhanced with as many additional unknowns as the number of edges along $\Gamma$, $N_{e}^{\Gamma}$.

To properly define this quadrature rule, we first obtain the inner product on each $T\in\mathcal{T}_{h}$ by mapping to the reference element $\hat{T}$. Using the Piola transformation, for any $\mathbf{q}$, $\mathbf{v}\in V_{h}$ and $\hat{\mathbf{q}}$, $\hat{\mathbf{v}}\in\hat{V}_{h}(\hat{T})$, we have

[TABLE]

If we choose $G|_{T}=J_{T}B_{T}^{-T}B_{T}^{-1}$ on each element $T$, we simplify the interaction of the functions on $\hat{T}$, thus obtaining $(G\mathbf{q},\mathbf{v})_{T}=(\hat{\mathbf{q}},\hat{\mathbf{v}})_{\hat{T}}$. Note that $G$ is indeed symmetric and positive definite for every $T\in\mathcal{T}_{h}$. In virtue of the previous result, the quadrature rule on $T$ is defined as (cf. [1, 2])

[TABLE]

where $|\hat{T}|$ is the area of $\hat{T}$ and $\hat{\mathbf{x}}_{c}=\frac{1}{3}\,(\hat{\mathbf{r}}_{1}+\hat{\mathbf{r}}_{2}+\hat{\mathbf{r}}_{3})$. This quadrature rule is exact for polynomials of degree $1$. Furthermore, if we denote by $\hat{\mathbf{v}}_{i}$ the basis function of $\hat{V}_{h}(\hat{T})$ associated with the $i$-th edge $\hat{e}_{i}$, for $i=1,2,3$, the following orthogonality condition is satisfied:

[TABLE]

where $|e_{i}|$ is the length of the corresponding $i$-th edge $e_{i}$ of $T$, for $i=1,2,3$. The global quadrature rule is thus given by

[TABLE]

Due to (6), the application of this quadrature rule permits to reduce $S$ to a diagonal matrix.

Summarizing, if we choose $\tilde{V}_{h}=V_{h}$ and further use the quadrature rule (7) for computing the integrals on the left-hand sides of (5b) and (5c), the semidiscrete scheme is now given by (5a), (5d) and (5e), in combination with the following equations:

[TABLE]

In [4], the unenhanced variant of this method (i.e., that obtained for a single subdomain) was shown to achieve $\mathcal{O}(h)$ optimal convergence for both pressures and velocities. In addition, the pressure variable was experimentally observed to be $\mathcal{O}(h^{2})$ superconvergent at the centroids of the triangles. Similar results were derived in [2] for the enhanced mixed method in the elliptic case. As reported in [10], on equilateral three-line meshes and if $K$ is the identity matrix, the normal velocity also achieves $\mathcal{O}(h^{2})$ superconvergence at the midpoints of the edges in the enhanced case. However, the order of convergence for the normal velocity is slightly less than 2 if a general three-line mesh and/or a general tensor $K$ are considered. The application of a post-processing technique further permits to increase the order of convergence of the velocity vector field to almost 2 (cf. [10]). In Section 5, all these results are confirmed through a collection of numerical experiments.

3.4 Matrix formulation

In this subsection, we describe how to express the expanded MFE scheme (5a), (8a), (8b), (5d) and (5e) in matrix form. Let us denote by $\{\mathbf{v}_{i}\}_{i=1}^{N_{V}}$, $\{w_{i}\}_{i=1}^{N_{W}}$ and $\{\mu_{i}\}_{i=1}^{N_{L}}$ the basis functions of $V_{h}$, $W_{h}$ and $L_{h}$, respectively. Then, the semidiscrete solution $(\mathbf{u}_{h},\tilde{\mathbf{u}}_{h},p_{h},\mathbf{\lambda}_{h})$ can be expressed as

[TABLE]

Omitting the time dependences, if we further define the vectors

[TABLE]

the differential system obtained above can be written as

[TABLE]

where the matrices $A_{1},A_{2}\in\mathbb{R}^{N_{V}\times N_{V}}$, $B\in\mathbb{R}^{N_{W}\times N_{V}}$ and $C\in\mathbb{R}^{N_{V}\times N_{L}}$ are given by

[TABLE]

By construction, $A_{1}$ and $A_{2}$ are block-diagonal matrices with $s$ diagonal blocks, each associated to a subdomain $\Omega_{i}$, for $i=1,2,\ldots,s$. In particular,

[TABLE]

where $A_{1}^{i},A_{2}^{i}\in\mathbb{R}^{(N_{V}/s)\times(N_{V}/s)}$. Note that the value $N_{V}/s\in\mathbb{N}$ corresponds to the number of edges per subdomain. Furthermore, $A_{1}$ is symmetric, positive definite and sparse, and $A_{2}$ is diagonal with positive diagonal entries (see (6)). In turn, $B$ is a rectangular block matrix of the form

[TABLE]

whose blocks $B^{i}\in\mathbb{R}^{(N_{W}/s)\times(N_{V}/s)}$, for $i=1,2,\ldots,s$. In this case, the value $N_{W}/s\in\mathbb{N}$ corresponds to the number of elements per subdomain. On the other hand, the diagonal matrix $D\in\mathbb{R}^{N_{W}\times N_{W}}$ is given by

[TABLE]

where $|T_{i}|$ denotes the area of $T_{i}$, for $i=1,2,\ldots,N_{W}$. Finally, the vectors $G_{h}^{D}\in\mathbb{R}^{N_{V}}$, $F_{h}\in\mathbb{R}^{N_{W}}$ and $G_{h}^{N}\in\mathbb{R}^{N_{L}}$ are defined to be

[TABLE]

The initial condition $P_{h}(0)=P_{h}^{0}\in\mathbb{R}^{N_{W}}$ is obtained as

[TABLE]

The preceding system (9) can be reduced to a system for ${P}_{h}$ and $\Lambda_{h}$ by eliminating both $U_{h}$ and $\tilde{U}_{h}$. To this end, we express

[TABLE]

where we introduce the submatrix

[TABLE]

and the vectors $\mathcal{B}=(0,B)$, $\mathcal{C}=(0,C)^{T}$, $\Phi_{h}=(\tilde{U}_{h},U_{h})^{T}$ and $R_{h}=(0,G_{h}^{D})^{T}$. From the first equation in (13), we obtain

[TABLE]

Inserting this expression into the second and third equations, the system can be written in the Schur complement form

[TABLE]

where the block elements of the system matrix are given by

[TABLE]

and the right-hand side terms are

[TABLE]

Note that $M\in\mathbb{R}^{N_{W}\times N_{W}}$ and $N\in\mathbb{R}^{N_{V}\times N_{V}}$ are both symmetric and positive definite. Moreover, $M$ is a block-diagonal matrix of the form

[TABLE]

where each block $M^{i}\in\mathbb{R}^{(N_{W}/s)\times(N_{W}/s)}$ is associated to a subdomain $\Omega_{i}$, for $i=1,2,\ldots,s$, and can be obtained as

[TABLE]

As a consequence, $M$ decouples across subdomains. Moreover, each block $M^{i}$ is a sparse matrix with, at most, 10 non-zero entries per row. The specific expressions for the coefficients of the 10-point stencil corresponding to $MP_{h}$ are derived in detail in Appendix A. Figures 11 and 12 further show the stencil of a pressure unknown associated to an upward and a downward triangle, respectively.

4 The fully discrete scheme

The resulting system (14) can be formulated in the form

[TABLE]

where

[TABLE]

This is a system of semi-explicit DAEs, namely: the system of ODEs (15a) subject to the constraint (15b). In this context, $P_{h}(t)$ is usually referred to as the differential component, whereas $\Lambda_{h}(t)$ is called the algebraic component. Since $\det(N)\neq 0$, (15) is indeed an index-1 DAE system. The initial condition for $P_{h}(t)$ is provided by (12), and that for $\Lambda_{h}(t)$ can be derived from (15b) as

[TABLE]

The initial values are thus consistent with (15), i.e., $F_{2}\left(0,P_{h}^{0},\Lambda_{h}^{0}\right)=0$.

It is well known that the time integration of index-1 DAE systems requires the application of stiffly accurate Runge–Kutta methods333A Runge–Kutta method is called stiffly accurate if its coefficients in the Butcher tableau satisfy $a_{qj}=b_{j}$, for $j=1,2,\ldots,q$. in order to obtain the same order of convergence for both the differential and algebraic components (cf. [17]). For instance, the families of Lobatto IIIA and Lobatto IIIC methods are convergent of order $2q-2$ for both components, $q$ being the number of internal stages of the corresponding method (cf. [16]). The simplest member in the Lobatto IIIA family (i.e., that corresponding to $q=2$) is the Crank–Nicolson method. Applied to (15), this method reads

[TABLE]

yielding approximations $P_{h}^{n}\approx P_{h}(t_{n})$ and $\Lambda_{h}^{n}\approx\Lambda_{h}(t_{n})$ on the equidistant time grid $0=t_{0}<t_{1}<\ldots<t_{n_{f}}=t_{f}$, where $t_{n}=n\tau$ and $\tau=t_{f}/n_{f}$, for $n=0,1,\ldots,n_{f}\in\mathbb{N}$. The notations $S_{h}^{n}=S_{h}(t_{n})$ and $T_{h}^{n}=T_{h}(t_{n})$ are also used.

In order to solve the system (16), we define $\hat{M}=I+\frac{1}{2}\,\tau D^{-1}M$ and consider the following procedure:

Set $P_{h}^{0}$ and $\Lambda_{h}^{0}$. 2. 2.

For $n=0,1,\ldots,n_{f}-1$:

1. (a)

Solve the system for $\Lambda_{h}^{n+1}$:

[TABLE] 2. (b)

Solve the system for $P_{h}^{n+1}$:

[TABLE]

The solution of (2a) may be obtained by using iterative linear solvers that require the computation of matrix-vector products involving the system matrix $N-\frac{1}{2}\,\tau Q^{T}\hat{M}^{-1}D^{-1}Q$. In such a case, we need to solve problems of the form

[TABLE]

for certain right-hand sides $b$. Taking into account that the coefficient matrix $\hat{M}$ decouples across subdomains, this equation involves independent subdomain problems which can be solved simultaneously. In addition, the use of hierarchical grids entails the same size for such subdomain problems, thus yielding a perfectly balanced workload among parallel processors. Once $\Lambda_{h}^{n+1}$ has been computed from (2a), the same ideas can be applied to the solution of (2b).

5 Numerical experiments

In this section, we study the numerical behaviour of the proposed method in the solution of a collection of parabolic initial-boundary value problems of the form (1): on the one hand, we examine its convergence properties by considering problems with a known analytical solution; on the other, we analyze its qualitative performance when applied to non-stationary flow models in porous media.

5.1 Convergence examples with known analytical solutions

5.1.1 A smooth solution test

Let us consider (1) posed on the irregular polygon with 7 sides shown in Figure 1. The polygon vertices are located at the points $(0,0)^{T}$, $(1,1)^{T}$, $(2.4,0.9)^{T}$, $(3.2,0)^{T}$, $(2.8,-0.8)^{T}$, $(1.9,-1.2)^{T}$ and $(0.7,-0.9)^{T}$. We further consider $t_{f}=2$, $\Gamma_{D}=\partial\Omega$, and

[TABLE]

The functions $f$, $p_{0}$ and $g_{D}$ are defined in such a way that

[TABLE]

is the exact solution of the problem. The spatial domain $\Omega$ is decomposed into $s=9$ subdomains, as shown in Figure 1 (left).

In the sequel, the pressure and velocity errors are computed by combining the $\ell^{\infty}$-norm in time with various norms in space. In particular, with an abuse of notation, the pressure errors are obtained in the norms

[TABLE]

and $\|r_{h}p-P_{h}\|_{\ell^{\infty}(\ell^{\infty})}$, where the maximum norm is used in both time and space. In these expressions, for any given $t$, $r_{h}p(t)\in\mathbb{R}^{N_{W}}$ is a vector whose $i$-th component is $p(\textbf{x}_{c}^{i},t)$, $\textbf{x}_{c}^{i}$ being the centroid of the $i$-th triangle of $\mathcal{T}_{h}$, for $i=1,2,\ldots,N_{W}$. Regarding the velocity variable, we consider the following norms

[TABLE]

In both expressions, $\mathbf{u}_{h}^{n}$ is a function of $V_{h}$ defined as $\mathbf{u}_{h}^{n}(\mathbf{x})=\sum_{i=1}^{N_{V}}U_{h,i}^{n}\,\mathbf{v}_{i}(\mathbf{x}),$ where $U_{h,i}^{n}$ are the elements of a vector $U_{h}^{n}$ containing the normal components of the numerical flux, given by

[TABLE]

In addition, the term $\left(\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n}\right)\cdot\mathbf{n}$ on the right-hand side of (21) denotes a vector in $\mathbb{R}^{N_{V}}$, whose elements are the normal components of $\mathbf{u}(t_{n})-\mathbf{u}_{h}^{n}$ at the midpoints of the edges of $\mathcal{T}_{h}$. The expression (20) further includes a linear operator $\mathbf{R}:V_{h}\rightarrow L^{2}(\Omega)$ that provides, for any $\mathbf{u}_{h}^{n}\in V_{h}$, a post-processed flux $\mathbf{R}(\mathbf{u}_{h}^{n})$ as defined in [10]. Specifically, $\mathbf{R}(\mathbf{u}_{h}^{n})$ is a piecewise function given by

[TABLE]

on each interior triangle $T\in\mathcal{T}_{h}$. The coefficients $a,b,\ldots,f$ are obtained as the solution of an overdetermined system with 9 equations by a least squares procedure. Such equations are obtained by imposing that the normal components of both the numerical and post-processed fluxes coincide at the midpoints of 9 edges: the three edges of $T$ and the two other edges of each of the three triangles bordering $T$. As described in [10], for those triangles containing an edge on the boundary of $\Omega$, a suitable modification of this procedure can be applied in order to preserve the global accuracy. In that work, the authors show theoretically that, for MFE methods, this post-processing technique is second-order accurate on three-line meshes if $K$ is the identity matrix. For enhanced MFE methods, numerical evidence reveals that the post-processed flux achieves close to second-order accuracy, even if a full tensor $K$ is considered.

The integral on the right-hand side of (20) is approximated element-wise by the midpoint quadrature rule. In turn, we use the formula

[TABLE]

for computing the integrals (11) and (12) involving the functions $f$ and $p_{0}$, respectively. Here, $\mathbf{x}_{m}^{i}$ denotes the midpoint of the $i$-th edge of $T$, for $i=1,2,3$. This formula is exact for polynomials of degree 2. Finally, the line integrals (10) for the Dirichlet boundary condition $g_{D}$ are approximated by Simpson’s rule.

Tables 1-4 show the pressure and velocity errors for a nested collection of spatial meshes $\mathcal{T}_{h}^{\ell}$, with $\ell=3,4,5,6$, and several decreasing time steps $\tau$. In particular, Tables 1 and 2 display the pressure errors computed with the norms $\|\cdot\|_{\ell^{\infty}(\ell^{2})}$ and $\|\cdot\|_{\ell^{\infty}(\ell^{\infty})}$, respectively, as defined above. In both cases, we can observe an unconditionally stable behaviour of the algorithm, which converges irrespective of the size of the parameters $h$ and $\tau$ under consideration. Furthermore, as revealed by the last row in both tables, the ratios of subsequent errors imply second-order convergence in space. Accordingly, the first errors in the last column show second-order convergence in time. Note that the last two positions in this column are meaningless, since they virtually represent the error due to the spatial discretization. On the other hand, Table 3 displays the computed errors for the post-processed fluxes according to formula (20). In this case, the method shows an unconditionally convergent behaviour, with second-order convergence in time and close to second-order convergence in space. Finally, in Table 4, we observe the errors for the normal fluxes given by (21). Once again, the method is unconditionally convergent; specifically, the order of convergence approaches an asymptotic value of 2 in time and approximately $1.6$ in space, in accordance with the numerical results reported in [10].

5.1.2 The case of discontinuous coefficients

As a second example, we consider a test problem inspired by a model by Mackinnon and Carey (cf. [18]) that involves a discontinuous permeability tensor. In particular, let (1) be posed on the unit square, with $t_{f}=3$ and $\Gamma_{D}=\partial\Omega$. In this case, $K(x,y)=k(x,y)\,I$, where $I$ is the $2\times 2$ identity matrix and $k(x,y)=k_{1}=1$, for $x\leq\frac{1}{2}$, and $k(x,y)=k_{2}=2$, for $x>\frac{1}{2}$. The data functions $f$, $p_{0}$ and $g_{D}$ are defined in such a way that the piecewise quadratic function

[TABLE]

is the exact solution of the problem, where

[TABLE]

The domain $\Omega$ is decomposed into $s=4$ triangular subdomains which are aligned with the permeability discontinuity, as shown in Figure 4. The integrals arising in the computations are approximated by the same quadrature rules as in the previous example.

Table 5 shows the pressure and velocity errors obtained for a nested collection of spatial meshes $\mathcal{T}_{h}^{\ell}$, with $\ell=1,2,3,4,5$, and a time step $\tau=1$E-01, considering the norms defined above. The use of the Crank–Nicolson time integrator for an exact solution of the form (22) permits us to obtain negligible errors in time and to study the spatial convergence of the method. The last row in the table contains the average orders of convergence. As expected, the method is second-order convergent for pressures and close to second-order convergent for post-processed and normal fluxes.

5.2 Numerical examples of time-dependent Darcy flow

5.2.1 A flow domain with low-permeability regions

In this example, we consider (1) posed on the unit square, with $t_{f}=5$, $f(x,y,t)=0$ and $p_{0}(x,y)=1-x$. Furthermore, $\Gamma_{D}=\{0,1\}\times(0,1)$ and $\Gamma_{N}=(0,1)\times\{0,1\}$. The pressure is specified to be equal to 1 on the left boundary and equal to 0 on the right boundary. A zero-flux condition is imposed on $\Gamma_{N}$. The flow domain contains two low-permeability regions, namely, $R_{1}=(0.2,0.3)\times(0,0.8)$ and $R_{2}=(0.6,0.7)\times(0.3,1)$. In particular, the permeability tensor is $K=k(x,y)\,I$, where $I$ is the $2\times 2$ identity matrix and

[TABLE]

This numerical example permits us to test the behaviour of the algorithm on problems involving abrupt variations in the permeability. A stationary version of this problem was considered in [12] in the framework of heat flow applications.

Figure 5 (left) shows the coarse grid $\mathcal{T}_{H}$ used in this example. It consists of $s=16$ triangular subdomains and it is adapted to the geometry of the low-permeability regions, depicted in the plot by shaded areas. In turn, Figure 5 (right) shows the fine grid $\mathcal{T}_{h}\equiv\mathcal{T}_{h}^{3}$ obtained after three regular refinement processes.

In Figure 6 (left), we display the pressure distribution obtained once the stationary state is reached. In this case, we consider the triangular mesh $\mathcal{T}_{h}^{5}$, obtained after five regular refinement processes of the coarse grid $\mathcal{T}_{H}$, and a time step $\tau=5$E-02. This figure illustrates the effect of the low-permeability regions on the pressure distribution. On the other hand, Figure 6 (right) shows the velocity field obtained at the stationary state. As usual, the length of the arrows is proportional to the module of the vectors. In this case, we consider the triangular mesh $\mathcal{T}_{h}^{3}$ shown in Figure 5 (right) and a time step $\tau=5$E-02. As expected from the physical configuration, no flow enters the low-permeability regions.

5.2.2 A flow domain with holes

Finally, we consider (1) posed on $\Omega\times(0,t_{f}]$, where $\Omega$ is a non-simply connected domain with two rectangular holes, and $t_{f}=5$. Specifically, we define $\Omega=(0,1)^{2}\setminus(R_{1}\cup R_{2})$, where $R_{1}=[0.2,0.6]\times[0.6,0.8]$ and $R_{2}=[0.4,0.8]\times[0.2,0.4]$. The source/sink term is $f(x,y,t)=0$, and the initial condition is $p_{0}(x,y)=1-x$. The boundary of $\Omega$ is decomposed into $\Gamma_{D}=\{0,1\}\times(0,1)$ and $\Gamma_{N}=(0,1)\times\{0,1\}\cup\partial R_{1}\cup\partial R_{2}$. In this context, the pressure is set to be equal to 1 on the boundary $\{0\}\times(0,1)$ and equal to 0 on the boundary $\{1\}\times(0,1)$. In turn, zero flux is imposed on $\Gamma_{N}$. In this example, we consider a uniform and isotropic permeability tensor $K=I$, where $I$ is the $2\times 2$ identity matrix. This experiment permits us to test the behaviour of the method on problems posed on non-simply connected domains. Stationary versions of this problem considering various configurations of holes in a square domain can be found in [12].

Figure 7 (left) shows the coarse grid $\mathcal{T}_{H}$ used in this example. It consists of $s=20$ triangular subdomains, which are aligned with the two rectangular holes. Figure 7 (right) shows the fine grid $\mathcal{T}_{h}\equiv\mathcal{T}_{h}^{3}$ obtained after three regular refinement processes.

The computed pressures and fluxes are displayed on Figure 8. The left-hand plot shows the pressure distribution at the stationary state, when considering the triangular mesh $\mathcal{T}_{h}^{5}$ and a time step $\tau=5$E-02. The effect of both holes on the pressure distribution is quite significant. In turn, the right-hand plot shows the velocity field obtained at the stationary state, when considering the triangular mesh $\mathcal{T}_{h}^{3}$ depicted in Figure 7 (right) and $\tau=5$E-02. In this case, the fluid flows around the internal boundaries of the domain.

6 Concluding remarks

A novel algorithm for solving time-dependent Darcy problems with possibly discontinuous tensor coefficients has been designed and tested on different geometries. The method is conceived to fit in the framework of hierarchical grids, taking advantage of their flexibility to handle complex problems while preserving the efficiency of stencil-based operations.

A relevant feature of the numerical procedure is that, provided that the material properties are constant within each subdomain, the corresponding stencil remains unchanged regardless of the level of refinement. In addition, the fully discrete scheme yields a local problem on each subdomain whose solution does not depend on the neighbouring subdomains, thus allowing for parallelization. Since the subdomain problems are all the same size, a perfectly balanced workload can be assigned among parallel processors in the solution process.

The structure of the resulting linear systems can lead to the development of very efficient linear solvers based on a judicious combination of preconditioning techniques and multigrid methods. Furthermore, the ease of implementation and robust numerical behaviour of the algorithm can make it suitable for more general flow problems. Both topics are the scope of current research.

Appendix A Derivation of the stencil for $MP_{h}$

In this appendix, we describe how to deduce the coefficients of the local 10-point stencil associated to $MP_{h}$, when $K$ is a tensor whose components do not depend on the spatial variables. Since matrix $M$ decouples across subdomains (cf. Subsection 3.4), it suffices to study each subdomain separately.

Let us consider a triangular subdomain $\Omega_{k}$, and let $\mathcal{T}_{h,k}\equiv\mathcal{T}_{h,k}^{\ell}$ be the conforming mesh that covers $\Omega_{k}$ with $4^{\ell}$ similar triangular elements. First of all, we need to set a suitable local numbering for the vertices, edges and triangles of the mesh. For a fixed refinement level $\ell$, each grid point in $\mathcal{T}_{h,k}$ is assigned a pair of indices $(i,j)$ that take the values $j=1,2,\ldots,2^{\ell}+1$ and $i=j,j+1,\ldots,2^{\ell}+1$. Based on this notation, the vertices of $\Omega_{k}$ are the points $(1,1)$, $(2^{\ell}+1,1)$ and $(2^{\ell}+1,2^{\ell}+1)$. In fact, such indices can be interpreted as the coordinates of the grid points in an oblique coordinate system $\{\mathbf{e}_{1}^{\prime},\mathbf{e}_{2}^{\prime}\}$ that depends on the geometry of the subdomain. Further, the edges associated to the $(i,j)$-th grid node are denoted as $e_{i,j}^{m}$, where $m=1,2,3$. This latter index determines which side of the coarse triangle $\Omega_{k}$ is parallel to the edge under consideration. Figure 9 shows the local numbering of a mesh $\mathcal{T}_{h,k}^{2}$ covering $\Omega_{k}$ that consists of 16 triangles. The three edges associated to the grid point $(3,2)$, namely $e_{3,2}^{1}$, $e_{3,2}^{2}$ and $e_{3,2}^{3}$, are depicted in blue, red and brown, respectively. Note that the triangles conforming the mesh can be classified, depending on its orientation, into upward and downward triangles. As shown in Figure 10, we will use the notation $T_{i,j}^{U}$ and $T_{i,j}^{D}$ to refer to the upward and downward triangles, respectively, that share the grid nodes $(i,j)$ and $(i+1,j+1)$. Figures 11 and 12 further show the distribution of several of these triangles on a certain region of a given mesh.

In the remaining of this section, we derive the stencil coefficients associated to an interior444In this framework, a triangle in $\mathcal{T}_{h,k}$ is said to be interior if the 10-point stencil of the pressure unknown associated to such a triangle is contained in the subdomain $\Omega_{k}$. upward triangle $T_{i,j}^{U}$. We will assume, in the sequel, that the unit vectors normal to the edges of $\mathcal{T}_{h,k}$ are defined to point outwards from $T_{i,j}^{U}$ and inwards to $T_{i,j}^{D}$.

To begin with, let us consider $w=w_{i,j}^{U}$ in the equation (5a), where $w_{i,j}^{U}$ is the basis function of $W_{h,k}$ associated to $T_{i,j}^{U}$ (i.e., the characteristic function of such a triangle). If we apply the divergence theorem on $T_{i,j}^{U}$ and take into account the orientation of normal unit vectors on $T_{i,j}^{U}$, we obtain

[TABLE]

where $P_{i,j}^{U}$ denotes the pressure unknown associated to the triangle $T_{i,j}^{U}$. In turn, $U_{i,j}^{2}$, $U_{i,j}^{3}$ and $U_{i+1,j}^{1}$ stand for the components of vector $U_{h}$ corresponding to the edges $e_{i,j}^{2}$, $e_{i,j}^{3}$ and $e_{i+1,j}^{1}$, respectively. As shown in Figure 10 (left), these are the sides of $T_{i,j}^{U}$. Their respective lengths are denoted by $l_{1}$, $l_{2}$ and $l_{3}$. Note that, due to the hierarchical structure of the grid covering $\Omega_{k}$, the edges $e_{i,j}^{m}$ all have the same length $l_{m}$, for $m=1,2,3$, regardless of the grid point coordinates $(i,j)$.

Next, we take $\mathbf{v}=\mathbf{v}_{i,j}^{2}$ in the equation (8b), where $\mathbf{v}_{i,j}^{2}$ is the basis function of $V_{h,k}$ associated to the edge $e_{i,j}^{2}$. Taking into account that $\mathbf{v}_{i,j}^{2}$ is different from zero only on the triangles $T_{i,j}^{U}$ and $T_{i,j}^{D}$ (cf. Figure 10 (left)), and using the quadrature rule introduced in (7), we get

[TABLE]

where $\tilde{U}_{i,j}^{m}$ denotes the component of vector $\tilde{U}_{h}$ associated to the edge $e_{i,j}^{m}$. Note that the five components of vector $\tilde{U}_{h}$ contained in the previous formula refer to the sides of the triangles $T_{i,j}^{U}$ and $T_{i,j}^{D}$, as shown in Figure 10 (left).

In order to compute the integrals on the right-hand side of (24), we first map each of them onto the reference element $\hat{T}$. Specifically, we introduce an affine mapping $F_{T}$, defined from the reference triangle $\hat{T}$, with vertices $\hat{\mathbf{r}}_{1}=(-1,0)^{T}$, $\hat{\mathbf{r}}_{2}=(1,0)^{T}$ and $\hat{\mathbf{r}}_{3}=(0,\sqrt{3})^{T}$, onto any triangle $T\in\mathcal{T}_{h,k}$, with vertices $\mathbf{r}_{1}=(x_{1},y_{1})^{T}$, $\mathbf{r}_{2}=(x_{2},y_{2})^{T}$ and $\mathbf{r}_{3}=(x_{3},y_{3})^{T}$, i.e.,

[TABLE]

where

[TABLE]

Then, we consider the first integral on the right-hand side of (24) and denote $T\equiv T_{i,j}^{U}$. Using the fact that $(\mathbf{v}_{i,j}^{2}\circ F_{T})(\hat{\mathbf{x}})=J_{T}^{-1}B_{T}\hat{\mathbf{v}}_{i,j}^{2}(\hat{\mathbf{x}})$ and $G|_{T}=J_{T}B_{T}^{-T}B_{T}^{-1}$, we obtain

[TABLE]

This same procedure can be applied to the rest of integrals in (24). In the sequel, we show that the integral on the right-hand side of the preceding expression is independent of the particular (upward or downward) triangle under consideration.

Let us consider the local ordering of vertices displayed on Figure 10 (right) for either an upward or a downward triangle. As the spatial mesh $\mathcal{T}_{h,k}$ is composed of similar triangles, it is easy to see that all mappings from $\hat{T}$ to an upward (or downward) triangle share the same Jacobian matrix $B_{T^{U}}$ (or $B_{T^{D}}$). Since $B_{T^{U}}=-B_{T^{D}}$, if tensor $K$ is assumed to be constant on $\Omega_{k}$, then $B_{T^{U}}^{-1}KB_{T^{U}}^{-T}=B_{T^{D}}^{-1}KB_{T^{D}}^{-T}$ for all the triangles in $\mathcal{T}_{h,k}$. We denote this product by $B_{k}^{-1}KB_{k}^{-T}$, and further define $J_{k}\equiv J_{T^{U}}=J_{T^{D}}$. In addition, it can be shown that $\hat{\mathbf{v}}_{i,j}^{m}=\frac{l_{m}}{2}\,\hat{\mathbf{v}}_{m}$, for any grid point $(i,j)$ and $m=1,2,3$, where

[TABLE]

are the standard basis functions of $\hat{V}_{h}(\hat{T})$. Therefore, introducing the notation555Note that, by construction, $a_{\ell,m}=a_{m,\ell}$.

[TABLE]

for $\ell,m\in\{1,2,3\}$, it is easy to derive, from (24), the following expression for $U_{i,j}^{2}$

[TABLE]

Accordingly, if we consider $\mathbf{v}=\mathbf{v}_{i,j}^{3}$ in the equation (8b), we will be able to express $U_{i,j}^{3}$ in terms of $\tilde{U}_{i,j}^{3}$, $\tilde{U}_{i,j-1}^{1}$, $\tilde{U}_{i+1,j}^{1}$, $\tilde{U}_{i,j}^{2}$ and $\tilde{U}_{i,j-1}^{2}$. Finally, if we take $\mathbf{v}=\mathbf{v}_{i+1,j}^{1}$ in (8b), $U_{i+1,j}^{1}$ will be given in terms of $\tilde{U}_{i+1,j}^{1}$, $\tilde{U}_{i,j}^{2}$, $\tilde{U}_{i+1,j}^{2}$, $\tilde{U}_{i,j}^{3}$ and $\tilde{U}_{i+1,j+1}^{3}$.

To conclude, we use the equation (8a) to express the components of vector $\tilde{U}_{h}$ as differences of pressure unknowns. In particular, let us consider $\mathbf{v}=\mathbf{u}_{i,j}^{2}$ in (8a). If we make use of the quadrature rule (7) and recall the orientation of normal unit vectors, we obtain

[TABLE]

In this expression, the component of $\tilde{U}_{h}$ corresponding to the edge $e_{i,j}^{2}$ is written in terms of the pressure unknowns associated to the triangles that share this edge. As shown in Figure 10 (left), such triangles are $T_{i,j}^{U}$ and $T_{i,j}^{D}$. If we repeat this procedure for the remaining 8 components of $\tilde{U}_{h}$ involved in the expressions above, we finally obtain 10 pressure unknowns: 7 of them are associated to upward triangles and the other 3 correspond to downward triangles (depicted in brown and blue, respectively, in Figure 11).

As a result, if we insert into (23) the expressions for $U_{i,j}^{m}$ in terms of $\tilde{U}_{i,j}^{m}$ derived above, for $m=1,2,3$ (see (26) for the case $U_{i,j}^{2}$), and further replace such $\tilde{U}_{i,j}^{m}$ by their corresponding formulas in terms of pressure differences (see (27) for the case $\tilde{U}_{i,j}^{2}$), we obtain an equation for the pressure variable at $T_{i,j}^{U}$ involving 10 unknowns, as shown in Figure 11. The term $MP_{h}$ arising in such an equation can be explicitly written using the standard notation for the elements of a stencil

[TABLE]

More precisely,

[TABLE]

where

[TABLE]

The entries $a_{\ell,m}$ in both matrices are given by (25). Remarkably, the stencil coefficients are the same for all interior pressures associated to upward triangles. Furthermore, these coefficients depend on tensor $K$ and the geometry of the given subdomain $\Omega_{k}$, but not on the number of mesh refinement processes.

The procedure to derive the stencil coefficients associated to an interior downward triangle is completely analogous. In such a case, we obtain

[TABLE]

where $S^{DD}=S^{UU}$, and

[TABLE]

Figure 12 shows the resulting stencil involving 10 pressure unknowns: 7 of them associated to downward triangles and the other 3 corresponding to upward triangles, coloured in blue and brown, respectively.

Note that, if we consider non-interior triangles, the coefficients for the pressure stencil can be derived following the same strategy; the only difference is that, in this case, the matrices $S^{UU}$, $S^{UD}$, $S^{DU}$ and $S^{DD}$ may contain a greater number of null elements.

Acknowledgments. This work was partially supported by MINECO grant MTM2014-52859.