- Email: [email protected]
A modified local and parallel finite element method for the mixed Stokes–Darcy model
J. Math. Anal. Appl. 435 (2016) 1129–1145
Contents lists available at ScienceDirect
Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa
A modified local and parallel finite element method for the mixed Stokes–Darcy model ✩ Guangzhi Du a , Yanren Hou a,b,∗ , Liyun Zuo a a b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Center for Computational Geosciences, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 17 May 2014 Available online 10 November 2015 Submitted by J. Guermond Keywords: Stokes equations Darcy’s law Decoupling Two-grid algorithm Parallel finite element method
a b s t r a c t Based on two-grid discretizations, a modified local and parallel finite element algorithm for the mixed Stokes–Darcy problem is proposed and analyzed in this paper. The idea is that, on the coarse grid, we solve a coupled Stokes–Darcy problem; on the fine grid, we first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve a series of independent local Stokes subproblems with the numerical solution of Darcy problem on fine grid to the interface coupling conditions in parallel. Optimal error estimates are obtained, and numerical experiments demonstrate the accuracy and efficiency of the modified local and parallel algorithm. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The coupling of incompressible fluid flow with porous media flow has received a wide publicity over the past years due to its wide applications in many practical problems and industrial settings, such as the interaction between surface and subsurface flows [12,13,17,24], the subsurface flow in karst aquifers [7,8,22], oil flow in vuggy porous [1], industrial filtrations [19], and so on. The incompressible fluid flow is described by the Stokes (or Navier–Stokes) equations, while the porous media flow is governed by the Darcy equations, two flows are coupled through certain interface conditions. A great deal of effort has been done to develop the effective methods of solving the mixed Stokes (or Navier–Stokes) and Darcy model, including domain decomposition methods [3,9,12,14,13,11,15,16,23], Lagrange multiplier methods [18,24], interface relaxation methods [27,29], precondition techniques [5,28], two-grid and multi-grid methods [4,6,10,30,39,38], partitioned time stepping algorithms [31,34,33], and so on. ✩ Subsidized by NSFC (Grant No. 11571274 & 11171269) and the Ph.D. Programs Foundation of Ministry of Education of China (Grant No. 20110201110027). * Corresponding author. E-mail addresses: [email protected] (G. Du), [email protected] (Y. Hou), [email protected] (L. Zuo).
http://dx.doi.org/10.1016/j.jmaa.2015.11.003 0022-247X/© 2015 Elsevier Inc. All rights reserved.
1130
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
The two-grid method was first successfully introduced to solve the mixed Stokes–Darcy problem by Mu and Xu in [30], then Cai and coworkers proposed a decoupled and linearized two-grid algorithm to solve a coupled Navier–Stokes/Darcy model in [6]; Chidyagwai and his/her coworker presented a two-grid scheme based on combinations of the continuous finite element method and the discontinuous Galerkin method in [10]. However, they all didn’t obtain the optimal error estimates for ∇u and p in L2 normal. In [39,38], a modified two-grid method is developed for the mixed Stokes–Darcy model, and the optimal error estimates are obtained. The idea of the modified two-grid method is that: on the coarse grid, solve a coupled Stokes–Darcy problem; on the fine grid, first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve the Stokes problem with the numerical solution of Darcy problem on fine grid to the interface coupling conditions. However, using the above two-grid algorithms, we find that on the fine grid the Stokes subproblem takes much more time than the Darcy problem. In order to further improve the effectiveness of solving the mixed Stokes–Darcy problem, we show a local and parallel finite element algorithm based on two-grid discretizations. It has been proposed for a class of linear and nonlinear elliptic boundary value problems by Xu and Zhou [35–37]. Then it was further extended to the Stokes equations [21], Navier–Stokes equations [20,26] and time-dependent convection–diffusion equations [25]. The idea of these algorithms is that the low frequency components can be approximated well by a relatively coarse grid and high frequency components can be computed on a fine grid by some local and parallel procedure. In this paper, following the idea of modified two-grid method in [39,38] and local and parallel approach in [21], we propose a modified local and parallel algorithm to solve the mixed Stokes–Darcy problem. The idea of our method is that: on the coarse grid, we solve a coupled Stokes–Darcy problem; on the fine grid, we first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve concurrently a series of independent local Stokes subproblems with the numerical solution of Darcy problem on fine grid to the interface coupling conditions. We obtain the optimal error estimates and numerical results show our scheme is efficient and effective. The rest of the paper is organized as follows. A mixed Stokes–Darcy model is described in the next section. In Section 3, a modified local and parallel algorithm is given. In Section 4, local a priori estimate is proved and convergence is deduced. Finally, some numerical experiments are presented in Section 5. 2. Model problem We consider the model in a bounded domain Ω ⊂ Rd (d = 2 or 3), which consists of a fluid region Ωf and a porous media region Ωp . Here Ωf ∩ Ωp = ∅, Ωf ∪ Ωp = Ω and ∂Ωf ∩ ∂Ωp = Γ is the interface. In addition, we denote Γf = ∂Ωf \ Γ, Γp = ∂Ωp \ Γ. The Stokes equations govern the fluid flow in Ωf : −νΔu + ∇p = f1
in Ωf ,
(2.1)
∇·u=0
in Ωf .
(2.2)
Here, u and p are the fluid velocity and pressure, ν is the kinetic viscosity, f1 is the external force. The Darcy equations govern the porous media flow in Ωp : up = −K∇ϕ
in Ωp ,
(2.3)
∇ · up = f2
in Ωp .
(2.4)
Here up and ϕ are the Darcy velocity and the piezometric head, K = {Kij }d×d is hydraulic conductivity tensor, denoting permeability of the rock, f2 is the source term. In this paper, we assume K is a symmetric and positive definite matrix.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1131
Combining Darcy’s law (2.3) with the continuity equation (2.4), we get the following equation: −∇ · (K∇ϕ) = f2
in Ωp .
(2.5)
On the interface Γ, we impose the following three interface conditions: u · nf + up · np = 0
on Γ,
(2.6)
∂u = gϕ on Γ, ∂nf ∂u νg u · τi 1 ≤ i ≤ (d − 1), =α −ντi ∂nf tr(K)
p − νnf
(2.7) on Γ,
(2.8)
where, nf and np denote the unit outward normal vectors on ∂Ωf and ∂Ωp , in particularly, nf = −np on the interface Γ. g is the gravitational acceleration, α is a positive parameter depending on the properties of the porous medium, τi , i = 1, . . . , d − 1, are the unit tangential vectors on the interface Γ. Eq. (2.8) is referred to as the Beavers–Joseph–Saffman interface condition [32], which is the simplified Beavers–Joseph interface condition [2]. For simplicity, we assume homogeneous Dirichlet boundary conditions are satisfied on Γf and Γp : u = 0 on Γf , ϕ=0
on Γp .
(2.9) (2.10)
The exact boundary conditions chosen above are not essential to either the analysis or the algorithms studied herein. For a domain E in Rd , the L2 (E) inner produce and norm for scalar and vector-valued functions are denoted by (·, ·)E and · 0,E , the norm of the Hilbert spaces H k (E) are denoted by · k,E . Define the following spaces: Hf = {v ∈ (H 1 (Ωf ))d : v = 0 on Γf }, Hp = {ψ ∈ H 1 (Ωp ) : ψ = 0 on Γp }, W = Hf × Hp , Q = L2 (Ωf ). The space W is equipped with the following norm: ∀ u = (u, ϕ) ∈ W , u0 = (u20,Ωf + ϕ20,Ωp )1/2 .
(2.11)
We shall assume that g, ν and α are constants, and f1 and f2 are smooth enough. Then the weak formulation of the above Stokes–Darcy model reads as follows: find u = (u, ϕ) ∈ W and p ∈ Q, such that
where
a(u, v) + b(v, p) = (f , v)
∀ v = (v, ψ) ∈ W,
(2.12)
b(u, q) = 0
∀ q ∈ Q,
(2.13)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1132
a(u, v) = af (u, v) + ap (ϕ, ψ) + aΓ (u, v), d−1 νg (u · τi )(v · τi ), af (u, v) = ν(∇u, ∇v)Ωf + α tr(K) i=1 Γ
ap (ϕ, ψ) = g(K∇ϕ, ∇ψ)Ωp , aΓ (u, v) = g (ϕv · nf − ψu · nf ), Γ
b(v, p) ≡ b(v, p) = −(p, ∇ · v)Ωf , (f , v) = (f1 , v)Ωf + g(f2 , ψ)Ωp . In [12] and [14], the well-posedness of the mixed Stokes–Darcy model (2.12)–(2.13) can be found. Throughout this paper, we shall use the letter c (with or without subscripts) to denote a generic positive constant which may depend the physical parameters and is independent of the mesh parameters. It may stand for different values at its different occurrences. Moreover, for the needs of theoretical analysis, we shall recall the following Poincaré and trace inequalities: there exist constants cp , ct that only depend on the domain Ωf , and cp , ct which only depend on the domain Ωp , such that for all v ∈ Hf and ψ ∈ Hp , v0,Ωf ≤ cp v1,Ωf ,
ψ0,Ωp ≤ cp ψ1,Ωp ,
(2.14)
vL2 (Γ) ≤ ct v1,Ωf ,
ψL2 (Γ) ≤ ct ψ1,Ωp .
(2.15)
3. Numerical algorithms We consider the regular triangulation τh (Ω) = {T } of Ω with mesh-size function h(x) whose value is the diameter hT of the element T containing x. We assume the triangulations τh (Ωf ) and τh (Ωp ) induced on the sub-domains Ωf and Ωp are compatible on the interface Γ. For the triangulations τh (Ωf ), a basic assumption [36] is that it is not exceedingly over-refined locally, i.e. A0. There exists γ ≥ 1 such that hγΩf ≤ ch(x)
∀x ∈ Ωf ,
(3.1)
where hΩf = max h(x) is the largest mesh size of τh (Ωf ). Sometimes, we will drop the subscript in hΩf and x∈Ωf
use h for the mesh size on a domain that is clear from the context. Moreover, we assume any subdomain of Ωf aligns with τh (Ωf ). Let Wh = Hf,h × Hp,h ⊂ W and Qh ⊂ Q be the finite element subspaces. The finite element spaces Hf,h and Qh approximating velocity and pressure in the fluid region are assumed to satisfy the discrete inf-sup condition: there exists a positive constant β, independent of h, such that ∀ vh ∈ Hf,h , qh ∈ Qh , b(vh , qh ) ≥ βvh 1,Ωf qh Q .
(3.2)
Given G ⊂⊂ Ω0 ⊂ Ωf (for G ⊂ Ω0 ⊂ Ωf , we write G ⊂⊂ Ω0 to mean that dist(∂G \∂Ωf , ∂Ω0 \∂Ωf ) > 0), and define Hf,h (G), Qh (G) and τh (G) to be the restriction of Hf,h , Qh and τh (Ωf ) to G, respectively, and 0 Hf,h (G) = {v ∈ Hf,h : supp v ⊂⊂ G}, Q0h (G) = {q ∈ Qh : supp q ⊂⊂ G}.
Similar to [21], we give some basic assumptions on the mixed finite element spaces in the fluid region.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1133
A1. Approximation. For each (u, p) ∈ H t+1 (G)d × H t (G)(t ≥ 1), then there exists an approximation (πh u, ρh p) ∈ Hf,h (G) × Qh (G) such that h−1 (u − πh u)0,G + u − πh u1,G ≤ chsG u1+s,G , 0 ≤ s ≤ t, h−1 (p − ρh p)−1,G + p − ρh p0,G ≤ chsG p1+s,G , 0 ≤ s ≤ t. A2. Inverse estimate. For any (v, q) ∈ Hf,h (G) × Qh (G), there holds v1,G ≤ ch−1 v0,G , q0,G ≤ ch−1 q−1,G . A3. Superapproximation. For G ⊂ Ω0 , let ω ∈ C0∞ (Ωf ) with supp ω ⊂⊂ G. Then for any (u, p) ∈ 0 Hf,h (G) × Qh (G), there is (v, q) ∈ Hf,h (G) × Q0h (G), such that h−1 (ωu − v)1,G ≤ cu1,G , h−1 (ωp − q)0,G ≤ cp0,G . There are several well-know finite element spaces satisfying the above assumptions, such as the MINI finite element and the P2 − P0 finite elements for t = 1, and the Hood–Taylor elements and the augmented P2 − P1 elements for t = 2. In this paper, we restrict our analysis to the case t = 1. Then the finite element scheme for (2.12)–(2.13) based on (Wh , Qh ) reads: Algorithm 1 (Coupling scheme): find uh = (uh , ϕh ) ∈ Wh and ph ∈ Qh such that a(uh , vh ) + b(vh , ph ) = (f , vh )
∀ vh = (vh , ψh ) ∈ Wh ,
(3.3)
b(uh , qh ) = 0
∀ qh ∈ Qh .
(3.4)
For the system (3.3)–(3.4), one has to solve a coupled problem, which usually result in difficulties in solving the resulting linear systems. In [30], Mu and Xu proposed a two-grid scheme to decouple the mixed 3 Stokes–Darcy problem, and the convergence rates of ∇u and p in L2 normal are O(h +H 2 ). Then a modified two-grid scheme in [39,38] is proposed and the optimal error estimates O(h + H 2 ) are obtained. Algorithm 2 (Modified two-grid scheme): Step 1. Solve the coupled problem on a coarse grid: find uH = (uH , ϕH ) ∈ WH ⊂ Wh and pH ∈ QH ⊂ Qh such that a(uH , vH ) + b(vH , pH ) = (f , vH )
∀ vH = (vH , ψH ) ∈ WH ,
(3.5)
b(uH , qH ) = 0
∀ qH ∈ QH .
(3.6)
Step 2. Solve a discrete Darcy problem in Ωp on the fine grid: find ϕh ∈ Hp,h such that ψh uH · nf ,
ap (ϕh , ψh ) = g(f2 , ψh )Ωp + g
∀ ψh ∈ Hp,h .
(3.7)
Γ
Step 3. On the fine grid, solve the discrete Stokes problem in Ωf : find uh ∈ Hf,h , ph ∈ Qh such that af (uh , vh ) + b(vh , ph ) = (f1 , vh )Ωf − g
ϕh vh · nf
∀ vh ∈ Hf,h ,
(3.8)
∀ qh ∈ Qh .
(3.9)
Γ h
b(u , qh ) = 0
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1134
We note that on the fine grid the Stokes problem takes much more time than the Darcy problem. To further improve the effectiveness, we combine the modified two-grid method and the local and parallel method [21] and propose a modified local and parallel scheme. Divide Ωf into a number of disjoint subdomains D1 , . . . , Dm , denote ΓDj = ∂Dj ∩ Γ = ∅. Then enlarge each Dj to Ωj such that Dj ⊂⊂ Ωj ⊂⊂ Ωf , and ΓΩj = ∂Ωj ∩ Γ. Algorithm 3 (Modified local and parallel scheme): Step 1. On a coarse grid, solve the coupled problem (3.5)–(3.6) to obtain uH = (uH , ϕH ). Step 2. On the fine grid, in the porous media region Ωp , solve the discrete Darcy problem (3.7) to get ϕh . Step 3. On the fine grid, in the fluid region Ωf , solve a series of independent local subproblems in parallel: find local fine grid corrections (ejh , ηhj ) ∈ Hf,h (Ωj ) × Qh (Ωj ) (j = 1, 2, . . . , m), for all (vh , qh ) ∈ Hf,h (Ωj ) × Qh (Ωj ) such that af (ejh , vh ) + b(vh , ηhj ) = (f1 , vh )Ωj − (af (uH , vh ) + b(vh , pH )) − g
ϕh vh · nf ,
(3.10)
ΓΩj
b(ejh , qh ) = −b(uH , qh ),
(3.11)
and set (uh , ph ) = (uH + ejh , pH + ηhj ) in Dj . 4. Theoretical analysis In this section, we shall first present local a priori error estimate for the Stokes finite element discretization on general shape regular grids. It will play a crucial role in our analysis. Then we consider the convergence of the modified local and parallel scheme we proposed. 4.1. Local a priori error estimate For a subdomain D of Ωf , it exists Ω0 such that D ⊂⊂ Ω0 ⊂ Ωf . To analyze the error estimate on a subdomain D, we need the following technical result. Lemma 4.1. Let ω ∈ C0∞ (Ωf ) such that supp ω ⊂⊂ Ω0 . Then ωw21,Ωf ≤ caf (w, ω 2 w) + cw20,Ω0 , ∀w ∈ Hf .
(4.1)
Proof. From the identity of af (·, ·) νωw21,Ωf ≤ af (ωw, ωw) = af (w, ω 2 w) + ν(∇ωw, ∇ωw)Ωf ≤ af (w, ω 2 w) + cw20,Ω0 , that is ωw21,Ωf ≤ caf (w, ω 2 w) + cw20,Ω0 .
2
We shall now present local a priori estimate for the Stokes finite element approximation that will be useful in our analysis. This type of estimates is an extension of the results in [21]. Lemma 4.2. Suppose that f ∈ L2 (ΓΩ0 ) and D ⊂⊂ Ω0 ⊂ Ωf . If Assumptions A0 and A1–A3 hold and (w, r) ∈ Hf,h × Qh satisfies
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1135
af (w, v) + b(v, r) =
f v,
∀v ∈ Hf,h (Ω0 ),
ΓΩ0
b(w, q) = 0,
∀q ∈ Qh (Ω0 ),
then w1,D + r0,D ≤ c(w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.2)
Proof. Let s be an integer such that s ≥ max{2γ, γ + 1}, Dj and Ωj (j = 1, 2, . . . , s) satisfy D1 ⊂⊂ D2 ⊂⊂ · · · ⊂⊂ Di ⊂⊂ · · · ⊂⊂ Ds ⊂⊂ Ωs , Ωs ⊂⊂ Ωs−1 ⊂⊂ · · · ⊂⊂ Ωj ⊂⊂ · · · ⊂⊂ Ω1 ⊂⊂ Ω0 . Choose G ⊂ Ωf satisfying D ⊂⊂ G ⊂⊂ D1 , let ω ∈ C0∞ (Ωf ) such that ω ≡ 1 on G and supp ω ⊂⊂ D1 . Note that r0,D ≤ ωr0,Ωf ≤ c
(∇ · φ, ωr) . φ1,Ωf
0 While for Ψ ∈ Hf,h (Ω0 ),
(∇ · φ, ωr) = (ω∇ · φ, r) = (∇ · (ωφ), r) − (φ∇ω, r)
= (∇ · (ωφ − Ψ), r) + af (w, Ψ) −
f Ψ − (φ∇ω, r).
ΓΩ0 0 According to (A3), choose Ψ ∈ Hf,h (D1 ), such that
ωφ − Ψ1,D1 ≤ chΩ0 φ1,D1 , then we have, Ψ1,D1 ≤ c(hΩ0 φ1,Ω0 + ωφ1,D1 ) ≤ cφ1,Ωf , |(∇ · (ωφ − Ψ), r)| ≤ cωφ − Ψ1,D1 r0,D1 ≤ chΩ0 r0,D1 φ1,Ωf , |
|af (w, Ψ)| ≤ cw1,D1 φ1,Ωf , f Ψ| ≤ cf L2 (ΓΩ0 ) ΨL2 (ΓΩ0 ) ≤ cf L2 (ΓΩ0 ) Ψ1,Ω0 ≤ cf L2 (ΓΩ0 ) φ1,Ωf ,
ΓΩ0
|(φ∇ω, r)| ≤ cφ1,D1 r−1,D1 ≤ cr−1,D1 φ1,Ωf . Hence |(∇ · φ, ωr)| ≤ c(hΩ0 φ1,Ωf r0,D1 + (w1,D1 + r−1,D1 + f L2 (ΓΩ0 ) )φ1,Ωf ). Then, applying the above estimates to (4.3), we have r0,D ≤ c(hΩ0 r0,D1 + w1,D1 + r−1,D1 + f L2 (ΓΩ0 ) ).
(4.3)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1136
Similarly, we obtain r0,Di−1 ≤ c(hΩ0 r0,Di + w1,Di + r−1,Di + f L2 (ΓΩ0 ) ), i = 1, 2, . . . , s − 1, where D0 = D. Thus r0,D ≤ c(hs−1 Ω0 r0,D s−1 + w1,D s−1 + r−1,D s−1 + f L2 (ΓΩ0 ) ) ≤ c(r−1,Ds−1 + w1,Ds−1 + f L2 (ΓΩ0 ) ) r0,D1 ≤ c(r−1,Ωs + w1,Ωs + f L2 (ΓΩ0 ) ).
(4.4)
0 On the other hand, there exists (v, q) ∈ Hf,h (D1 ) × Q0h (D1 ) such that
ω 2 w − v1,D1 ≤ chΩ0 w1,D1 , ω 2 r − q0,D1 ≤ chΩ0 r0,D1 , which imply |af (w, ω 2 w − v)| ≤ chΩ0 w21,D1 ,
(4.5)
|(∇ · (v − ω 2 w), r) + (∇ · w, ω 2 r − q)| ≤ chΩ0 w1,D1 r0,D1 , f v| ≤ cf L2 (ΓD1 ) v1,D1 ≤ cf L2 (ΓΩ0 ) (hΩ0 w1,D1 + ωw1,Ωf ). |
(4.6) (4.7)
ΓD 1
Note that (∇ · (ω 2 w), r) = (∇ · w, ω 2 r) + (2ωw∇ω, r) = (∇ · w, ω 2 r − q) + (2ωw∇ω, r).
(4.8)
We have from Lemma 3.1 and (4.5)–(4.8) that c−1 ωw21,Ωf ≤ af (w, ω 2 w) + w20,Ω0
= af (w, ω w − v) + (∇ · v, r) + 2
f v + w20,Ω0
ΓD 1
= af (w, ω w − v) + (∇ · (ω w), r) + (∇ · (v − ω w), r) + 2
2
f v + w20,Ω0
2
ΓD 1
≤ c(hΩ0 (w21,D1 + w1,D1 r0,D1 ) + w20,Ω0 ) + c(ωw1,D1 r−1,D1 + f L2 (ΓΩ0 ) (hΩ0 w1,D1 + ωw1,Ωf )). Hence 1/2
ωw1,Ωf ≤ c(hΩ0 (w1,D1 + r0,D1 ) + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.9)
Using (4.4) and (4.9), we get 1/2
w1,D ≤ c(hΩ0 w1,Ωs + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.10)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1137
Similarly, 1/2
w1,Ωj ≤ c(hΩ0 w1,Ωj−1 + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ), j = 1, 2, . . . , s. Therefore s
w1,Ωs ≤ c(hΩ2 0 w1,Ω0 + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ) ≤ c(w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ). Combining (4.11) with (4.4), (4.10), we obtain (4.2).
(4.11)
2
4.2. Convergence analysis In this subsection, we consider the convergence of the modified local and parallel scheme. Moreover, we assume the solution (u, p, ϕ) of (2.1)–(2.8) satisfies u ∈ (H 2 (Ωf ))d , p ∈ H 1 (Ωf ) and ϕ ∈ H 2 (Ωp ). In addition, we define the following norms: |u − u |1,Ωf h
m 1 =( u − uh 21,Dj ) 2 ,
|p − p |0,Ωf h
j=1
m 1 =( p − ph 20,Dj ) 2 .
(4.12)
j=1
The following lemma about the convergence of the coupling scheme can be found in [30]. Lemma 4.3. For the problem (3.3)–(3.4), we have the following error estimate u − uh 1,Ωf + p − ph 0,Ωf + ϕ − ϕh 1,Ωp ≤ ch.
(4.13)
Denote e ≡ (u − uh , ϕ − ϕh ). To obtain e0 + p − ph −1,Ωf , we consider the dual problem of the Stokes–Darcy model: given g ∈ (L2 (Ωf ))d × L2 (Ωp ), θ ∈ L2 (Ωf ) ∩ H 1 (Ωf ), find (w, r) = (w, ζ, r) ∈ Hf × Hp × Q, such that ∀(v, q) = (v, ψ, q) ∈ Hf × Hp × Q, a(v, w) − b(v, r) + b(w, q) = (g, v) + (θ, q)Ωf .
(4.14)
Assume the solution of the problem (4.14) satisfies (w, ζ, r) ∈ (H 2 (Ωf ))d × H 2 (Ωp ) × H 1 (Ωf ), then we have the regularity 1
(w22,Ωf + ζ22,Ωp ) 2 + r1,Ωf ≤ c(g0 + θ1,Ωf ).
(4.15)
Moreover, there are the following approximation properties: there exists (wh , ζh , rh ) ∈ Hf,h × Hp,h × Qh such that w − wh 1,Ωf ≤ chw2,Ωf ,
(4.16)
ζ − ζh 1,Ωp ≤ chζ2,Ωp ,
(4.17)
r − rh 0,Ωf ≤ chr1,Ωf .
(4.18)
Lemma 4.4. Under the assumption of regularity (4.15), we have the following estimate 1
(u − uh 20,Ωf + ϕ − ϕh 20,Ωp ) 2 + p − ph −1,Ωf ≤ ch2 .
(4.19)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1138
Proof. Comparing the problems (2.12)–(2.13) and (3.3)–(3.4), then taking (vh , qh ) = (wh , rh ) = (wh , ζh , rh ) gives a(e, wh ) + b(wh , p − ph ) = 0,
(4.20)
b(u − uh , rh ) = 0.
(4.21)
Setting v = e, q = p − ph in (4.14), we have (g, e) + (θ, p − ph )Ωf = a(e, w) − b(u − uh , r) + b(w, p − ph ) = a(e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph ).
(4.22)
The right side in (4.22) are bounded by Schwarz and trace inequalities, |a(e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph )| = |af (u − uh , w − wh ) + ap (ϕ − ϕh , ζ − ζh ) + aΓ (e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph )| ≤ c(u − uh 1,Ωf (w − wh 1,Ωf + ζ − ζh 1,Ωp + r − rh 0,Ωf ) + ϕ − ϕh 1,Ωp (w − wh 1,Ωf + ζ − ζh 1,Ωp ) + p − ph 0,Ωf w − wh 1,Ωf ).
(4.23)
Combining (4.13), (4.16)–(4.18), (4.23) with (4.22) yields (g, e) + (θ, p − ph )Ωf ≤ ch2 (w2,Ωf + ζ2,Ωp + r1,Ωf ) 1
≤ ch2 ((w22,Ωf + ζ22,Ωp ) 2 + r1,Ωf ) ≤ ch2 (g0 + θ1,Ωf ),
(4.24)
(4.19) thus follows. 2 Next, we give the convergence of the modified two-grid scheme, it can be found in [39,38]. Lemma 4.5. Let (uh , ph , ϕh ) be the solution of modified two-grid scheme, then the following error estimates hold: u − uh 1,Ωf ≤ c(h + H 2 ),
(4.25)
p − ph 0,Ωf ≤ c(h + H 2 ),
(4.26)
ϕ − ϕh 1,Ωp ≤ c(h + H 2 ).
(4.27)
Finally, we will consider the convergence of the modified local and parallel scheme. Theorem 4.1. Let (uh , ph , ϕh ) be the solution of the modified local and parallel scheme, we have the following error estimates |u − uh |1,Ωf + |p − ph |0,Ωf ≤ c(h + H 2 ).
(4.28)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1139
Proof. To prove (4.28), we first prove: u − uh 1,Dj + p − ph 0,Dj ≤ c(h + H 2 ).
(4.29)
Let vh = (vh , 0) in (3.3), we have af (uh , vh ) + b(vh , ph ) = (f1 , vh )Ωf − g
ϕh vh · nf ,
∀vh ∈ Hf,h (Ωf )
(4.30)
∀qh ∈ Qh (Ωf ).
(4.31)
Γ
b(uh , qh ) = 0, Comparing the Algorithm 3 and (4.30)–(4.31), gives af (uh − uh , vh ) + b(vh , ph − ph ) = g
(ϕh − ϕh )vh · nf
∀vh ∈ Hf,h (Ωj ),
(4.32)
∀qh ∈ Qh (Ωj ).
(4.33)
ΓΩj
b(uh − uh , qh ) = 0 Note that from [30,39,38], we have ϕh − ϕh 1,Ωf ≤ cH 2 .
(4.34)
By Lemma 4.2, trace inequality (2.15), (4.34) and (4.19), we have uh − uh 1,Dj + ph − ph 0,Dj ≤ c(ϕh − ϕh L2 (ΓΩj ) + uh − uh 0,Ωj + ph − ph −1,Ωj ) ≤ c(ϕh − ϕh 1,Ωf + uh − uH 0,Ωj + ph − pH −1,Ωj + ejh 0,Ωj + ηhj −1,Ωj ) ≤ c(H 2 + ejh 0,Ωj + ηhj −1,Ωj ).
(4.35)
To estimate ejh 0,Ωj + ηhj −1,Ωj , we consider the following problem ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ∂w ⎪ −ντi ∂n ⎪ f ⎪ ⎪ ⎩
−νΔw − ∇r = θ in Ωj , ∇ · w = φ in Ωj , = g(ϕh − ϕH )w/v on ΓΩj ,
∂w r − νnf ∂n f
νg = α tr(K) w · τi 1 ≤ i ≤ (d − 1) on ΓΩj ,
(4.36)
w = 0 on ∂Ωj \ ΓΩj .
There holds w2,Ωj + r1,Ωj ≤ c(θ0,Ωj + φ1,Ωj ).
(4.37)
Then the weak form of (4.36) can be read as follows: for ∀(v, q) ∈ Hf (Ωj ) × L2 (Ωj ) af (w, v) − b(v, r) + b(w, q) + g ΓΩj
(ϕh − ϕH )w · nf = (θ, v)Ωj + (φ, q)Ωj .
(4.38)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1140
Let (wμ , rμ ) ∈ Hf,μ (Ωj ) × Qμ (Ωj ) be the finite element approximation of (w, r), i.e., af (w − wμ , v) − b(v, r − rμ ) + g
(ϕh − ϕH )(w − wμ ) · nf = 0 ∀v ∈ Hf,μ (Ωj ),
(4.39)
ΓΩj
b(w − wμ , q) = 0 ∀q ∈ Qμ (Ωj ),
(4.40)
w − wμ 1,Ωj + r − rμ 0,Ωj ≤ cμ(w2,Ωj + r1,Ωj ) ≤ cμ(θ0,Ωj + φ1,Ωj ).
(4.41)
where μ = h or H. Then
It follows from the problem (4.30)–(4.31) and (3.10)–(3.11) that af (ejh , v) + b(v, ηhj ) − b(ejh , q) = af (uh − uH , v) + b(v, ph − pH ) − b(uh − uH , q) +g (ϕh − ϕh )v · nf , ∀(v, q) ∈ Hf,h (Ωj ) × Qh (Ωj ), ΓΩj
(4.42)
af (uh − uH , wH ) + b(wH , ph − pH ) − b(uh − uH , rH ) + g
(ϕh − ϕH )wH · nf = 0.
(4.43)
ΓΩj
Setting v = ejh , q = ηhj in (4.38), then combining (4.39)–(4.40) with (4.42)–(4.43) yields (θ, ejh ) + (φ, ηhj ) =
af (ejh , w)
−
b(ejh , r)
+
b(w, ηhj )
(ϕh − ϕH )w · nf
+g ΓΩj
=
af (ejh , wh )
−
b(ejh , rh )
+
b(wh , ηhj )
+g
(ϕh − ϕH )wh · nf
ΓΩj
= af (uh − uH , wh ) + b(wh , ph − pH ) − b(uh − uH , rh ) + g
(2ϕh − ϕh − ϕH )wh · nf
ΓΩj
= af (uh − uH , wh − wH ) + b(wh − wH , ph − pH ) − b(uh − uH , rh − rH ) +g (ϕh − ϕH )(wh − wH ) · nf + g (ϕh − ϕh )wh · nf . ΓΩj
(4.44)
ΓΩj
Hence, we can derive from (4.41), (4.44), trace inequality (2.15), (4.13) and (4.34) that |(θ, ejh ) + (φ, ηhj )| ≤ c[H(uh − uH 1,Ωf + ph − pH 0,Ωf + ϕh − ϕH 1,Ωf ) + ϕh − ϕh 1,Ωf ](θ0,Ωj + φ1,Ωj ) ≤ cH 2 (θ0,Ωj + φ1,Ωj ),
(4.45)
which implies ejh 0,Ωj + ηhj −1,Ωj ≤ cH 2 .
(4.46)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1141
Applying (4.46) to (4.35), we get uh − uh 1,Dj + ph − ph 0,Dj ≤ cH 2 .
(4.47)
Then (4.29) follows from triangle inequalities and (4.13). Finally, from the definition of |u − uh |1,Ωf and |p − ph |0,Ωf in (4.12), there holds |u − uh |1,Ωf + |p − ph |0,Ωf ≤
m (u − uh 1,Dj + p − ph 0,Dj ) j=1
≤ c(h + H 2 ).
2
(4.48)
5. Numerical experiments In this section, we give three numerical experiments. In the first two experiments, we consider the cases for which the exact solution is known and unknown respectively, to demonstrate the accuracy and efficiency of the modified local and parallel method. In the final case, we choose different K to examine the effects of the hydraulic conductivity tensor. Assume the computational domain Ωf = (0, 1) × (1, 2) and Ωp = (0, 1) × (0, 1) with interface Γ = (0, 1)×{1}. To apply the modified local and parallel algorithm, we divide the domain Ωf into four subdomains D1 = (0, 1/4) × (1, 2),
D2 = (1/4, 1/2) × (1, 2),
D3 = (1/2, 3/4) × (1, 2),
D4 = (3/4, 1) × (1, 2),
and set Ω1 = (0, 5/16) × (1, 2),
Ω2 = (3/16, 9/16) × (1, 2),
Ω3 = (7/16, 13/16) × (1, 2),
Ω4 = (11/16, 1) × (1, 2).
The well-known Mini elements (P 1b − P 1) are used for the Stokes problem and the linear Lagrangian elements (P 1) for the Darcy flow. The relation of the coarse mesh H and the fine mesh h is chosen as h = H 2 . Moreover, we apply the software package FreeFEM++ to implement the algorithms and all results are generated by the same computer. For simplicity of notation, we denote ehϑ = ϑh − ϑ,
h
ehϑ = ϑ − ϑ,
where ϑ can be ϕ, u, p. Experiment 1. We set the physical parameters g, ν, α equal to 1, and K = I. The exact solution is given by ⎧ 2 2 2 3 ⎪ ⎨ u = ([x (y − 1) + y], [− 3 x(y − 1) + 2 − π sin(πx)]), p = [2 − π sin(πx)] sin( π2 y), ⎪ ⎩ ϕ = [2 − π sin(πx)][1 − y − cos(πy)]. The forcing terms and the boundary conditions follow the exact solution.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1142
Table 1 The convergence performance and CPU time of the coupling scheme. h
∇(ϕ − ϕh )0,Ωp
Rate
∇(u − uh )0,Ωf
Rate
p − ph 0,Ωf
Rate
CPU
1 16 1 36 1 64
0.375149 0.163785 0.092593
– 1.02 0.99
0.360524 0.162884 0.091482
– 0.98 1.00
0.148571 0.055708 0.033513
– 1.21 0.88
5.893 34.106 243.437
Table 2 The convergence performance and CPU time of the modified two-grid method. h
h ∇eϕ 0,Ωp
Rate
∇eh u 0,Ωf
Rate
eh p 0,Ωf
Rate
Stokes
Darcy
CPU
1 16 1 36 1 64 1 100
0.375473 0.166299 0.094427 0.060576
– 1.00 0.98 1.02
0.370676 0.161416 0.091266 0.058505
– 1.03 0.99 1.02
0.166929 0.056210 0.032502 0.019653
– 1.34 0.95 1.15
0.125 0.655 2.262 5.975
0.015 0.078 0.265 0.655
0.281 1.279 3.978 9.922
Table 3 The convergence performance and CPU time of the modified local and parallel scheme. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
Stokes
Darcy
CPU
1 16 1 36 1 64 1 100
0.375473 0.166299 0.094427 0.060576
− 1.00 0.98 1.02
0.370912 0.163275 0.091858 0.059705
− 1.01 1.00 0.99
0.159621 0.060694 0.029189 0.017481
− 1.19 1.27 1.18
0.062 0.25 0.78 2.044
0.015 0.093 0.266 0.655
0.232 0.89 2.511 6.038
In Table 1, we list the errors between the exact solution and the finite element approximate solutions of the coupled problem, convergent rates and computational time with the mesh size h = 1/16, 1/36, 1/64. It’s obvious that the errors for ∇ϕ, ∇u and p in the L2 -norm are all around of the order of O(h). In Table 2, the errors and convergence rates of the modified two-grid method are given. Observe that, the numerical results exhibit the first-order reduction. Meanwhile, the computational time of the Stokes problem (3.8)–(3.9), the Darcy problem (3.7), and the Stokes–Darcy problem (3.5)–(3.9) are listed. It is showed that the Stokes problem takes more time than the Darcy problem. In Table 3, we show the convergence performance and computational time of the modified local and parallel scheme with the mesh size h = 1/16, 1/36, 1/64, 1/100. Especially, “Stokes” denotes the max CPU time spent by each local Stokes sub-problem (3.10)–(3.11), and “Darcy” means the computational time of the problem (3.7). Apparently, the numerical results agree with the theoretical analysis strongly. That is to say, the modified local and parallel scheme has the same accuracy as the coupling scheme. From Table 2 and Table 3, it’s obvious that the modified local and parallel scheme spends less time than the modified two-grid scheme. Experiment 2. In this numerical experiment, we consider a case for which the exact solution for the coupled problem is unknown. The physical parameters are the same as in Experiment 1, the boundary conditions are homogeneous Dirichlet boundary conditions: u=0
on Γf ,
ϕ = 0 on Γp .
The external forces are f1 = (−(A1 + B1 ) + C1 , −(A2 + B2 ) + C2 ), f2 = −12(2x − 1)y, where
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1143
Table 4 The convergence performance of the modified local and parallel scheme. h
∇(ϕh0 − ϕh )0,Ωp
Rate
|∇(uh0 − uh )|0,Ωf
Rate
|ph0 − ph |0,Ωf
Rate
1 16 1 25 1 36 1 49 1 64
0.076086 0.047841 0.035011 0.029221 0.020406
– 1.04 1.01 0.96 1.00
0.070361 0.046777 0.034825 0.027308 0.020989
– 0.95 0.92 0.91 0.92
0.022657 0.017673 0.012593 0.010220 0.007295
– 0.99 0.98 0.94 0.97
Table 5 The convergence performance of the modified local and parallel algorithm with k = 2. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
0.017839 0.007818 0.004467 0.002838
– 1.02 0.97 1.02
0.045087 0.020268 0.011283 0.007215
– 0.99 1.02 1.00
0.016630 0.006532 0.003287 0.001909
– 1.15 1.19 1.22
Table 6 The convergence performance of the modified local and parallel algorithm with k = 0.2. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
0.178335 0.078160 0.044599 0.028389
– 1.02 0.98 1.01
0.045089 0.020269 0.011284 0.007215
– 0.99 1.02 1.00
0.021258 0.009331 0.004455 0.002741
– 1.02 1.28 1.09
A1 = 2(6x2 − 6x + 1)(2y 3 − 6y 2 + 5y − 2), B1 = 12x2 (x − 1)2 (y − 1), C1 = 2(6x2 − 6x + 1)y, A2 = −6(2x − 1)(y − 2)2 (y 2 + 1), B2 = −2x(x − 1)(2x − 1)(6y 2 − 12y + 5), C2 = 2x(x − 1)(2x − 1). In this case, we first solve the coupled Stokes–Darcy model in a fine mesh grid h0 (h0 = 1/64), derive the approximate solution (uh0 , ph0 , ϕh0 ), and regard the approximate solution as one “exact solution”. In Table 4, we show the errors between the “exact solution” and the modified local and parallel scheme solutions and give the convergence rate. Experiment 3. The final case is carried out in order to assess the influence of the physical parameters on the convergence rate of the modified local and parallel scheme. In this case, we take K = kI and choose the exact solution ⎧ 2 ⎪ ⎨ u = ((y − 1) , x(x − 1)), p = 2ν(x + y − 1), ⎪ ⎩ ϕ = k−1 (x(1 − x)(y − 1) + (y − 1)3 /3) + 2ν x. g The forcing terms and the boundary conditions follow the exact solution, the physical parameters g, ν, α are the same as in Experiment 1. In Tables 5–7, we list the convergence performance of the modified local and parallel algorithm with K = 2I, K = 0.2I, K = 0.02I, respectively. The errors of ∇ϕ and p in L2 -normal increase as k decreases,
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1144
Table 7 The convergence performance of the modified local and parallel algorithm with k = 0.02. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
1.783670 0.781297 0.446450 0.284395
– 1.02 0.97 1.01
0.045188 0.020327 0.011315 0.007235
– 0.99 1.02 1.00
0.112322 0.056715 0.030252 0.019952
– 0.84 1.09 0.93
while the errors of ∇u in L2 -normal almost remain the same. The changes of k do not influence the convergence rate of ∇ϕ, p and ∇u. On the whole, we proposed a modified local and parallel algorithm for the mixed Stokes and Darcy problem. Optimal error estimates are derived. Numerical experiments show that the modified local and parallel scheme we proposed not only has the same accuracy as the coupling scheme, but also saves a lot of time than the modified two-grid scheme. We thus conclude that the modified local and parallel algorithm is very effective. References [1] T. Arbogast, D.S. Brunson, A computational method for approximating a Darcy–Stokes system governing a vuggy porous medium, Comput. Geosci. 11 (2007) 207–218. [2] G. Beavers, D. Joseph, Boundary conditions at a naturally impermeable, J. Fluid Mech. 30 (1967) 197–207. [3] Y. Boubendir, S. Tlupova, Domain decomposition methods for solving Stokes–Darcy problems with boundary integrals, SIAM J. Sci. Comput. 35 (2013) B82–B106. [4] M.C. Cai, M. Mu, A multilevel decoupled method for a mixed Stokes/Darcy model, J. Comput. Appl. Math. 236 (2012) 2452–2465. [5] M.C. Cai, M. Mu, J.C. Xu, Preconditioning techniques for a mixed Stokes/Darcy model in porous media applications, J. Comput. Appl. Math. 233 (2009) 346–355. [6] M.C. Cai, M. Mu, J.C. Xu, Numerical solution to a mixed Navier–Stokes/Darcy model by the two-grid approach, SIAM J. Numer. Anal. 47 (2009) 3325–3338. [7] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, W. Zhao, Finite element approximation for Stokes–Darcy flow with Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 47 (2010) 4239–4256. [8] Y. Cao, M. Gunzburger, F. Hua, X. Wang, Coupled Stokes–Darcy model with Beavers–Joseph interface boundary condition, Commun. Math. Sci. 8 (2010) 1–25. [9] W. Chen, M. Gunzburger, F. Hua, X.M. Wang, A parallel Robin–Robin domain decomposition method for the Stokes– Darcy system, SIAM J. Numer. Anal. 49 (2011) 1064–1084. [10] P. Chidyagwai, B. Rivière, A two-grid method for coupled free flow with porous media flow, Adv. Water Resour. 34 (2011) 1113–1123. [11] M. Discaaaiati, A. Quarteroni, A. Valli, Robin–Robin domain decomposition methods for the Stokes–Darcy coupling, SIAM J. Numer. Anal. 45 (2007) 1246–1268. [12] M. Discacciati, Domain decomposition methods for the coupling of surface and groundwater flows, Ph.D. dissertation, École Polytechnique Fédérale de Lausanne, 2004. [13] M. Discacciati, A. Quarteroni, Convergence analysis of a subdomain iterative method for the finite element approximation of the coupling of Stokes and Darcy equations, Comput. Vis. Sci. 6 (2004) 93–103. [14] M. Discacciati, E. Miglio, A. Quarteroni, Mathematical and numerical models for coupling surface and groundwater flows, Appl. Numer. Math. 43 (2002) 57–74. [15] M. Discacciati, A. Quarteroni, Analysis of a domain decomposition method for the coupling Stokes and Darcy equations, in: F. Brezzi, et al. (Eds.), Numerical Analysis and Advanced Applications, Enumath 2001, Springer, Milan, 2003, pp. 3–20. [16] J. Galvis, M. Sarkis, Balancing domain decomposition methods for mortar coupling Stokes–Darcy Systems, in: Proceedings of the 16th International Conference on Domain Decomposition Methods, New York, 2005. [17] V. Girault, B. Rivière, DG approximation of coupled Navier–Stokes and Darcy equations by Beaver–Joseph–Saffman interface condition, SIAM J. Numer. Anal. 47 (2009) 2052–2089. [18] R. Glowinski, T. Pan, J. Periaux, A Lagrange multiplier/fictitious domain method for the numerical simulation of incompressible viscous flow around moving grid bodies: I. Case where the rigid body motions are known a priori, C. R. Acad. Sci. Paris Ser. I Math. 324 (1997) 361–369. [19] N. Hanspal, A. Waghode, V. Nassehi, R. Wakeman, Numerical analysis of coupled Stokes/Darcy flow in industrial filtrations, Transp. Porous Media 64 (2006) 73–101. [20] Y.N. He, J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms for the Navier–Stokes problem, J. Comput. Math. 24 (2006) 227–238. [21] Y.N. He, J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms for the Stokes problem, Numer. Math. 109 (2008) 415–434. [22] F. Hua, Modeling, analysis and simulation of Stokes–Darcy system with Beavers–Joseph interface condition, Ph.D. dissertation, The Florida State University, 2009.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1145
[23] B. Jiang, A parallel domain decomposition method for coupling of surface and groundwarter flows, Comput. Methods Appl. Mech. Engrg. 198 (2009) 947–957. [24] W.J. Layton, F. Schieweck, I. Yotov, Coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 40 (2003) 2195–2218. [25] Q.F. Liu, Y.R. Hou, Local and parallel finite element algorithms for time-dependent convection–diffusion equations, Appl. Math. Mech. (English Ed.) 30 (2009) 787–794. [26] F. Ma, Y. Ma, W. Wo, Local and parallel finite element algorithms based on two-grid discretization for steady Navier–Stokes equations, Appl. Math. Mech. (English Ed.) 28 (2007) 27–35. [27] S. Markus, E. Houstis, A. Catlin, J. Rice, P. Tsompanopoulou, E. Vavalis, D. Gottfried, K. Su, G. Balakrishnan, An agent-based netcentric framework for multidisciplinary problem solving environments (MPSE), Int. J. Comput. Eng. Sci. 1 (2000) 33–60. [28] A. Marquez, S. Meddahi, F.J. Sayas, A decoupled preconditioning technique for a mixed Stokes–Darcy model, J. Sci. Comput. 57 (2013) 174–192. [29] M. Mu, Solving composite problems with interface relaxation, SIAM J. Sci. Comput. 20 (1999) 1394–1416. [30] M. Mu, J.C. Xu, A two-grid method of a mixed Stokes–Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 45 (2007) 1801–1813. [31] M. Mu, X.H. Zhu, Decoupled schemes for a non-stationary mixed Stokes–Darcy model, Math. Comp. 79 (2010) 707–731. [32] P. Saffman, On the boundary condition at the surface of a porous media, Stud. Appl. Math. 50 (1971) 93–101. [33] L. Shan, H.B. Zheng, Partitioned time stepping method for fully evolutionary Stokes–Darcy flow with the Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 51 (2013) 813–839. [34] L. Shan, H.B. Zheng, W.J. Layton, A decoupling method with different subdomain time steps for the nonstationary Stokes–Darcy model, Numer. Methods Partial Differential Equations 29 (2013) 549–583. [35] J.C. Xu, A.H. Zhou, Some local and parallel properties of finite element discretizations, in: Proceedings for Eleventh International Conference on Domain Decomposition Methods, Greenwich, 1999, pp. 140–147. [36] J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms based on two-grid discretizations, Math. Comp. 69 (1999) 881–909. [37] J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms based on two-grid discretizations for nonlinear problems, Adv. Comput. Math. 14 (2001) 293–327. [38] T. Zhang, J.Y. Yuan, Two novel decoupling algorithms for the steady Stokes–Darcy model based on two-grid discretizations, Discrete Contin. Dyn. Syst. Ser. B 19 (2014) 849–865. [39] L.Y. Zuo, Y.R. Hou, A decoupling two-grid algorithm for the mixed Stokes–Darcy model with the Beavers–Joseph interface condition, Numer. Methods Partial Differential Equations 3 (2014) 1066–1082.
Contents lists available at ScienceDirect
Journal of Mathematical Analysis and Applications www.elsevier.com/locate/jmaa
A modified local and parallel finite element method for the mixed Stokes–Darcy model ✩ Guangzhi Du a , Yanren Hou a,b,∗ , Liyun Zuo a a b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China Center for Computational Geosciences, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
a r t i c l e
i n f o
Article history: Received 17 May 2014 Available online 10 November 2015 Submitted by J. Guermond Keywords: Stokes equations Darcy’s law Decoupling Two-grid algorithm Parallel finite element method
a b s t r a c t Based on two-grid discretizations, a modified local and parallel finite element algorithm for the mixed Stokes–Darcy problem is proposed and analyzed in this paper. The idea is that, on the coarse grid, we solve a coupled Stokes–Darcy problem; on the fine grid, we first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve a series of independent local Stokes subproblems with the numerical solution of Darcy problem on fine grid to the interface coupling conditions in parallel. Optimal error estimates are obtained, and numerical experiments demonstrate the accuracy and efficiency of the modified local and parallel algorithm. © 2015 Elsevier Inc. All rights reserved.
1. Introduction The coupling of incompressible fluid flow with porous media flow has received a wide publicity over the past years due to its wide applications in many practical problems and industrial settings, such as the interaction between surface and subsurface flows [12,13,17,24], the subsurface flow in karst aquifers [7,8,22], oil flow in vuggy porous [1], industrial filtrations [19], and so on. The incompressible fluid flow is described by the Stokes (or Navier–Stokes) equations, while the porous media flow is governed by the Darcy equations, two flows are coupled through certain interface conditions. A great deal of effort has been done to develop the effective methods of solving the mixed Stokes (or Navier–Stokes) and Darcy model, including domain decomposition methods [3,9,12,14,13,11,15,16,23], Lagrange multiplier methods [18,24], interface relaxation methods [27,29], precondition techniques [5,28], two-grid and multi-grid methods [4,6,10,30,39,38], partitioned time stepping algorithms [31,34,33], and so on. ✩ Subsidized by NSFC (Grant No. 11571274 & 11171269) and the Ph.D. Programs Foundation of Ministry of Education of China (Grant No. 20110201110027). * Corresponding author. E-mail addresses: [email protected] (G. Du), [email protected] (Y. Hou), [email protected] (L. Zuo).
http://dx.doi.org/10.1016/j.jmaa.2015.11.003 0022-247X/© 2015 Elsevier Inc. All rights reserved.
1130
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
The two-grid method was first successfully introduced to solve the mixed Stokes–Darcy problem by Mu and Xu in [30], then Cai and coworkers proposed a decoupled and linearized two-grid algorithm to solve a coupled Navier–Stokes/Darcy model in [6]; Chidyagwai and his/her coworker presented a two-grid scheme based on combinations of the continuous finite element method and the discontinuous Galerkin method in [10]. However, they all didn’t obtain the optimal error estimates for ∇u and p in L2 normal. In [39,38], a modified two-grid method is developed for the mixed Stokes–Darcy model, and the optimal error estimates are obtained. The idea of the modified two-grid method is that: on the coarse grid, solve a coupled Stokes–Darcy problem; on the fine grid, first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve the Stokes problem with the numerical solution of Darcy problem on fine grid to the interface coupling conditions. However, using the above two-grid algorithms, we find that on the fine grid the Stokes subproblem takes much more time than the Darcy problem. In order to further improve the effectiveness of solving the mixed Stokes–Darcy problem, we show a local and parallel finite element algorithm based on two-grid discretizations. It has been proposed for a class of linear and nonlinear elliptic boundary value problems by Xu and Zhou [35–37]. Then it was further extended to the Stokes equations [21], Navier–Stokes equations [20,26] and time-dependent convection–diffusion equations [25]. The idea of these algorithms is that the low frequency components can be approximated well by a relatively coarse grid and high frequency components can be computed on a fine grid by some local and parallel procedure. In this paper, following the idea of modified two-grid method in [39,38] and local and parallel approach in [21], we propose a modified local and parallel algorithm to solve the mixed Stokes–Darcy problem. The idea of our method is that: on the coarse grid, we solve a coupled Stokes–Darcy problem; on the fine grid, we first consider the discrete Darcy problem with a coarse grid approximation uH to the interface coupling conditions, then solve concurrently a series of independent local Stokes subproblems with the numerical solution of Darcy problem on fine grid to the interface coupling conditions. We obtain the optimal error estimates and numerical results show our scheme is efficient and effective. The rest of the paper is organized as follows. A mixed Stokes–Darcy model is described in the next section. In Section 3, a modified local and parallel algorithm is given. In Section 4, local a priori estimate is proved and convergence is deduced. Finally, some numerical experiments are presented in Section 5. 2. Model problem We consider the model in a bounded domain Ω ⊂ Rd (d = 2 or 3), which consists of a fluid region Ωf and a porous media region Ωp . Here Ωf ∩ Ωp = ∅, Ωf ∪ Ωp = Ω and ∂Ωf ∩ ∂Ωp = Γ is the interface. In addition, we denote Γf = ∂Ωf \ Γ, Γp = ∂Ωp \ Γ. The Stokes equations govern the fluid flow in Ωf : −νΔu + ∇p = f1
in Ωf ,
(2.1)
∇·u=0
in Ωf .
(2.2)
Here, u and p are the fluid velocity and pressure, ν is the kinetic viscosity, f1 is the external force. The Darcy equations govern the porous media flow in Ωp : up = −K∇ϕ
in Ωp ,
(2.3)
∇ · up = f2
in Ωp .
(2.4)
Here up and ϕ are the Darcy velocity and the piezometric head, K = {Kij }d×d is hydraulic conductivity tensor, denoting permeability of the rock, f2 is the source term. In this paper, we assume K is a symmetric and positive definite matrix.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1131
Combining Darcy’s law (2.3) with the continuity equation (2.4), we get the following equation: −∇ · (K∇ϕ) = f2
in Ωp .
(2.5)
On the interface Γ, we impose the following three interface conditions: u · nf + up · np = 0
on Γ,
(2.6)
∂u = gϕ on Γ, ∂nf ∂u νg u · τi 1 ≤ i ≤ (d − 1), =α −ντi ∂nf tr(K)
p − νnf
(2.7) on Γ,
(2.8)
where, nf and np denote the unit outward normal vectors on ∂Ωf and ∂Ωp , in particularly, nf = −np on the interface Γ. g is the gravitational acceleration, α is a positive parameter depending on the properties of the porous medium, τi , i = 1, . . . , d − 1, are the unit tangential vectors on the interface Γ. Eq. (2.8) is referred to as the Beavers–Joseph–Saffman interface condition [32], which is the simplified Beavers–Joseph interface condition [2]. For simplicity, we assume homogeneous Dirichlet boundary conditions are satisfied on Γf and Γp : u = 0 on Γf , ϕ=0
on Γp .
(2.9) (2.10)
The exact boundary conditions chosen above are not essential to either the analysis or the algorithms studied herein. For a domain E in Rd , the L2 (E) inner produce and norm for scalar and vector-valued functions are denoted by (·, ·)E and · 0,E , the norm of the Hilbert spaces H k (E) are denoted by · k,E . Define the following spaces: Hf = {v ∈ (H 1 (Ωf ))d : v = 0 on Γf }, Hp = {ψ ∈ H 1 (Ωp ) : ψ = 0 on Γp }, W = Hf × Hp , Q = L2 (Ωf ). The space W is equipped with the following norm: ∀ u = (u, ϕ) ∈ W , u0 = (u20,Ωf + ϕ20,Ωp )1/2 .
(2.11)
We shall assume that g, ν and α are constants, and f1 and f2 are smooth enough. Then the weak formulation of the above Stokes–Darcy model reads as follows: find u = (u, ϕ) ∈ W and p ∈ Q, such that
where
a(u, v) + b(v, p) = (f , v)
∀ v = (v, ψ) ∈ W,
(2.12)
b(u, q) = 0
∀ q ∈ Q,
(2.13)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1132
a(u, v) = af (u, v) + ap (ϕ, ψ) + aΓ (u, v), d−1 νg (u · τi )(v · τi ), af (u, v) = ν(∇u, ∇v)Ωf + α tr(K) i=1 Γ
ap (ϕ, ψ) = g(K∇ϕ, ∇ψ)Ωp , aΓ (u, v) = g (ϕv · nf − ψu · nf ), Γ
b(v, p) ≡ b(v, p) = −(p, ∇ · v)Ωf , (f , v) = (f1 , v)Ωf + g(f2 , ψ)Ωp . In [12] and [14], the well-posedness of the mixed Stokes–Darcy model (2.12)–(2.13) can be found. Throughout this paper, we shall use the letter c (with or without subscripts) to denote a generic positive constant which may depend the physical parameters and is independent of the mesh parameters. It may stand for different values at its different occurrences. Moreover, for the needs of theoretical analysis, we shall recall the following Poincaré and trace inequalities: there exist constants cp , ct that only depend on the domain Ωf , and cp , ct which only depend on the domain Ωp , such that for all v ∈ Hf and ψ ∈ Hp , v0,Ωf ≤ cp v1,Ωf ,
ψ0,Ωp ≤ cp ψ1,Ωp ,
(2.14)
vL2 (Γ) ≤ ct v1,Ωf ,
ψL2 (Γ) ≤ ct ψ1,Ωp .
(2.15)
3. Numerical algorithms We consider the regular triangulation τh (Ω) = {T } of Ω with mesh-size function h(x) whose value is the diameter hT of the element T containing x. We assume the triangulations τh (Ωf ) and τh (Ωp ) induced on the sub-domains Ωf and Ωp are compatible on the interface Γ. For the triangulations τh (Ωf ), a basic assumption [36] is that it is not exceedingly over-refined locally, i.e. A0. There exists γ ≥ 1 such that hγΩf ≤ ch(x)
∀x ∈ Ωf ,
(3.1)
where hΩf = max h(x) is the largest mesh size of τh (Ωf ). Sometimes, we will drop the subscript in hΩf and x∈Ωf
use h for the mesh size on a domain that is clear from the context. Moreover, we assume any subdomain of Ωf aligns with τh (Ωf ). Let Wh = Hf,h × Hp,h ⊂ W and Qh ⊂ Q be the finite element subspaces. The finite element spaces Hf,h and Qh approximating velocity and pressure in the fluid region are assumed to satisfy the discrete inf-sup condition: there exists a positive constant β, independent of h, such that ∀ vh ∈ Hf,h , qh ∈ Qh , b(vh , qh ) ≥ βvh 1,Ωf qh Q .
(3.2)
Given G ⊂⊂ Ω0 ⊂ Ωf (for G ⊂ Ω0 ⊂ Ωf , we write G ⊂⊂ Ω0 to mean that dist(∂G \∂Ωf , ∂Ω0 \∂Ωf ) > 0), and define Hf,h (G), Qh (G) and τh (G) to be the restriction of Hf,h , Qh and τh (Ωf ) to G, respectively, and 0 Hf,h (G) = {v ∈ Hf,h : supp v ⊂⊂ G}, Q0h (G) = {q ∈ Qh : supp q ⊂⊂ G}.
Similar to [21], we give some basic assumptions on the mixed finite element spaces in the fluid region.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1133
A1. Approximation. For each (u, p) ∈ H t+1 (G)d × H t (G)(t ≥ 1), then there exists an approximation (πh u, ρh p) ∈ Hf,h (G) × Qh (G) such that h−1 (u − πh u)0,G + u − πh u1,G ≤ chsG u1+s,G , 0 ≤ s ≤ t, h−1 (p − ρh p)−1,G + p − ρh p0,G ≤ chsG p1+s,G , 0 ≤ s ≤ t. A2. Inverse estimate. For any (v, q) ∈ Hf,h (G) × Qh (G), there holds v1,G ≤ ch−1 v0,G , q0,G ≤ ch−1 q−1,G . A3. Superapproximation. For G ⊂ Ω0 , let ω ∈ C0∞ (Ωf ) with supp ω ⊂⊂ G. Then for any (u, p) ∈ 0 Hf,h (G) × Qh (G), there is (v, q) ∈ Hf,h (G) × Q0h (G), such that h−1 (ωu − v)1,G ≤ cu1,G , h−1 (ωp − q)0,G ≤ cp0,G . There are several well-know finite element spaces satisfying the above assumptions, such as the MINI finite element and the P2 − P0 finite elements for t = 1, and the Hood–Taylor elements and the augmented P2 − P1 elements for t = 2. In this paper, we restrict our analysis to the case t = 1. Then the finite element scheme for (2.12)–(2.13) based on (Wh , Qh ) reads: Algorithm 1 (Coupling scheme): find uh = (uh , ϕh ) ∈ Wh and ph ∈ Qh such that a(uh , vh ) + b(vh , ph ) = (f , vh )
∀ vh = (vh , ψh ) ∈ Wh ,
(3.3)
b(uh , qh ) = 0
∀ qh ∈ Qh .
(3.4)
For the system (3.3)–(3.4), one has to solve a coupled problem, which usually result in difficulties in solving the resulting linear systems. In [30], Mu and Xu proposed a two-grid scheme to decouple the mixed 3 Stokes–Darcy problem, and the convergence rates of ∇u and p in L2 normal are O(h +H 2 ). Then a modified two-grid scheme in [39,38] is proposed and the optimal error estimates O(h + H 2 ) are obtained. Algorithm 2 (Modified two-grid scheme): Step 1. Solve the coupled problem on a coarse grid: find uH = (uH , ϕH ) ∈ WH ⊂ Wh and pH ∈ QH ⊂ Qh such that a(uH , vH ) + b(vH , pH ) = (f , vH )
∀ vH = (vH , ψH ) ∈ WH ,
(3.5)
b(uH , qH ) = 0
∀ qH ∈ QH .
(3.6)
Step 2. Solve a discrete Darcy problem in Ωp on the fine grid: find ϕh ∈ Hp,h such that ψh uH · nf ,
ap (ϕh , ψh ) = g(f2 , ψh )Ωp + g
∀ ψh ∈ Hp,h .
(3.7)
Γ
Step 3. On the fine grid, solve the discrete Stokes problem in Ωf : find uh ∈ Hf,h , ph ∈ Qh such that af (uh , vh ) + b(vh , ph ) = (f1 , vh )Ωf − g
ϕh vh · nf
∀ vh ∈ Hf,h ,
(3.8)
∀ qh ∈ Qh .
(3.9)
Γ h
b(u , qh ) = 0
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1134
We note that on the fine grid the Stokes problem takes much more time than the Darcy problem. To further improve the effectiveness, we combine the modified two-grid method and the local and parallel method [21] and propose a modified local and parallel scheme. Divide Ωf into a number of disjoint subdomains D1 , . . . , Dm , denote ΓDj = ∂Dj ∩ Γ = ∅. Then enlarge each Dj to Ωj such that Dj ⊂⊂ Ωj ⊂⊂ Ωf , and ΓΩj = ∂Ωj ∩ Γ. Algorithm 3 (Modified local and parallel scheme): Step 1. On a coarse grid, solve the coupled problem (3.5)–(3.6) to obtain uH = (uH , ϕH ). Step 2. On the fine grid, in the porous media region Ωp , solve the discrete Darcy problem (3.7) to get ϕh . Step 3. On the fine grid, in the fluid region Ωf , solve a series of independent local subproblems in parallel: find local fine grid corrections (ejh , ηhj ) ∈ Hf,h (Ωj ) × Qh (Ωj ) (j = 1, 2, . . . , m), for all (vh , qh ) ∈ Hf,h (Ωj ) × Qh (Ωj ) such that af (ejh , vh ) + b(vh , ηhj ) = (f1 , vh )Ωj − (af (uH , vh ) + b(vh , pH )) − g
ϕh vh · nf ,
(3.10)
ΓΩj
b(ejh , qh ) = −b(uH , qh ),
(3.11)
and set (uh , ph ) = (uH + ejh , pH + ηhj ) in Dj . 4. Theoretical analysis In this section, we shall first present local a priori error estimate for the Stokes finite element discretization on general shape regular grids. It will play a crucial role in our analysis. Then we consider the convergence of the modified local and parallel scheme we proposed. 4.1. Local a priori error estimate For a subdomain D of Ωf , it exists Ω0 such that D ⊂⊂ Ω0 ⊂ Ωf . To analyze the error estimate on a subdomain D, we need the following technical result. Lemma 4.1. Let ω ∈ C0∞ (Ωf ) such that supp ω ⊂⊂ Ω0 . Then ωw21,Ωf ≤ caf (w, ω 2 w) + cw20,Ω0 , ∀w ∈ Hf .
(4.1)
Proof. From the identity of af (·, ·) νωw21,Ωf ≤ af (ωw, ωw) = af (w, ω 2 w) + ν(∇ωw, ∇ωw)Ωf ≤ af (w, ω 2 w) + cw20,Ω0 , that is ωw21,Ωf ≤ caf (w, ω 2 w) + cw20,Ω0 .
2
We shall now present local a priori estimate for the Stokes finite element approximation that will be useful in our analysis. This type of estimates is an extension of the results in [21]. Lemma 4.2. Suppose that f ∈ L2 (ΓΩ0 ) and D ⊂⊂ Ω0 ⊂ Ωf . If Assumptions A0 and A1–A3 hold and (w, r) ∈ Hf,h × Qh satisfies
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1135
af (w, v) + b(v, r) =
f v,
∀v ∈ Hf,h (Ω0 ),
ΓΩ0
b(w, q) = 0,
∀q ∈ Qh (Ω0 ),
then w1,D + r0,D ≤ c(w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.2)
Proof. Let s be an integer such that s ≥ max{2γ, γ + 1}, Dj and Ωj (j = 1, 2, . . . , s) satisfy D1 ⊂⊂ D2 ⊂⊂ · · · ⊂⊂ Di ⊂⊂ · · · ⊂⊂ Ds ⊂⊂ Ωs , Ωs ⊂⊂ Ωs−1 ⊂⊂ · · · ⊂⊂ Ωj ⊂⊂ · · · ⊂⊂ Ω1 ⊂⊂ Ω0 . Choose G ⊂ Ωf satisfying D ⊂⊂ G ⊂⊂ D1 , let ω ∈ C0∞ (Ωf ) such that ω ≡ 1 on G and supp ω ⊂⊂ D1 . Note that r0,D ≤ ωr0,Ωf ≤ c
(∇ · φ, ωr) . φ1,Ωf
0 While for Ψ ∈ Hf,h (Ω0 ),
(∇ · φ, ωr) = (ω∇ · φ, r) = (∇ · (ωφ), r) − (φ∇ω, r)
= (∇ · (ωφ − Ψ), r) + af (w, Ψ) −
f Ψ − (φ∇ω, r).
ΓΩ0 0 According to (A3), choose Ψ ∈ Hf,h (D1 ), such that
ωφ − Ψ1,D1 ≤ chΩ0 φ1,D1 , then we have, Ψ1,D1 ≤ c(hΩ0 φ1,Ω0 + ωφ1,D1 ) ≤ cφ1,Ωf , |(∇ · (ωφ − Ψ), r)| ≤ cωφ − Ψ1,D1 r0,D1 ≤ chΩ0 r0,D1 φ1,Ωf , |
|af (w, Ψ)| ≤ cw1,D1 φ1,Ωf , f Ψ| ≤ cf L2 (ΓΩ0 ) ΨL2 (ΓΩ0 ) ≤ cf L2 (ΓΩ0 ) Ψ1,Ω0 ≤ cf L2 (ΓΩ0 ) φ1,Ωf ,
ΓΩ0
|(φ∇ω, r)| ≤ cφ1,D1 r−1,D1 ≤ cr−1,D1 φ1,Ωf . Hence |(∇ · φ, ωr)| ≤ c(hΩ0 φ1,Ωf r0,D1 + (w1,D1 + r−1,D1 + f L2 (ΓΩ0 ) )φ1,Ωf ). Then, applying the above estimates to (4.3), we have r0,D ≤ c(hΩ0 r0,D1 + w1,D1 + r−1,D1 + f L2 (ΓΩ0 ) ).
(4.3)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1136
Similarly, we obtain r0,Di−1 ≤ c(hΩ0 r0,Di + w1,Di + r−1,Di + f L2 (ΓΩ0 ) ), i = 1, 2, . . . , s − 1, where D0 = D. Thus r0,D ≤ c(hs−1 Ω0 r0,D s−1 + w1,D s−1 + r−1,D s−1 + f L2 (ΓΩ0 ) ) ≤ c(r−1,Ds−1 + w1,Ds−1 + f L2 (ΓΩ0 ) ) r0,D1 ≤ c(r−1,Ωs + w1,Ωs + f L2 (ΓΩ0 ) ).
(4.4)
0 On the other hand, there exists (v, q) ∈ Hf,h (D1 ) × Q0h (D1 ) such that
ω 2 w − v1,D1 ≤ chΩ0 w1,D1 , ω 2 r − q0,D1 ≤ chΩ0 r0,D1 , which imply |af (w, ω 2 w − v)| ≤ chΩ0 w21,D1 ,
(4.5)
|(∇ · (v − ω 2 w), r) + (∇ · w, ω 2 r − q)| ≤ chΩ0 w1,D1 r0,D1 , f v| ≤ cf L2 (ΓD1 ) v1,D1 ≤ cf L2 (ΓΩ0 ) (hΩ0 w1,D1 + ωw1,Ωf ). |
(4.6) (4.7)
ΓD 1
Note that (∇ · (ω 2 w), r) = (∇ · w, ω 2 r) + (2ωw∇ω, r) = (∇ · w, ω 2 r − q) + (2ωw∇ω, r).
(4.8)
We have from Lemma 3.1 and (4.5)–(4.8) that c−1 ωw21,Ωf ≤ af (w, ω 2 w) + w20,Ω0
= af (w, ω w − v) + (∇ · v, r) + 2
f v + w20,Ω0
ΓD 1
= af (w, ω w − v) + (∇ · (ω w), r) + (∇ · (v − ω w), r) + 2
2
f v + w20,Ω0
2
ΓD 1
≤ c(hΩ0 (w21,D1 + w1,D1 r0,D1 ) + w20,Ω0 ) + c(ωw1,D1 r−1,D1 + f L2 (ΓΩ0 ) (hΩ0 w1,D1 + ωw1,Ωf )). Hence 1/2
ωw1,Ωf ≤ c(hΩ0 (w1,D1 + r0,D1 ) + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.9)
Using (4.4) and (4.9), we get 1/2
w1,D ≤ c(hΩ0 w1,Ωs + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ).
(4.10)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1137
Similarly, 1/2
w1,Ωj ≤ c(hΩ0 w1,Ωj−1 + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ), j = 1, 2, . . . , s. Therefore s
w1,Ωs ≤ c(hΩ2 0 w1,Ω0 + w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ) ≤ c(w0,Ω0 + r−1,Ω0 + f L2 (ΓΩ0 ) ). Combining (4.11) with (4.4), (4.10), we obtain (4.2).
(4.11)
2
4.2. Convergence analysis In this subsection, we consider the convergence of the modified local and parallel scheme. Moreover, we assume the solution (u, p, ϕ) of (2.1)–(2.8) satisfies u ∈ (H 2 (Ωf ))d , p ∈ H 1 (Ωf ) and ϕ ∈ H 2 (Ωp ). In addition, we define the following norms: |u − u |1,Ωf h
m 1 =( u − uh 21,Dj ) 2 ,
|p − p |0,Ωf h
j=1
m 1 =( p − ph 20,Dj ) 2 .
(4.12)
j=1
The following lemma about the convergence of the coupling scheme can be found in [30]. Lemma 4.3. For the problem (3.3)–(3.4), we have the following error estimate u − uh 1,Ωf + p − ph 0,Ωf + ϕ − ϕh 1,Ωp ≤ ch.
(4.13)
Denote e ≡ (u − uh , ϕ − ϕh ). To obtain e0 + p − ph −1,Ωf , we consider the dual problem of the Stokes–Darcy model: given g ∈ (L2 (Ωf ))d × L2 (Ωp ), θ ∈ L2 (Ωf ) ∩ H 1 (Ωf ), find (w, r) = (w, ζ, r) ∈ Hf × Hp × Q, such that ∀(v, q) = (v, ψ, q) ∈ Hf × Hp × Q, a(v, w) − b(v, r) + b(w, q) = (g, v) + (θ, q)Ωf .
(4.14)
Assume the solution of the problem (4.14) satisfies (w, ζ, r) ∈ (H 2 (Ωf ))d × H 2 (Ωp ) × H 1 (Ωf ), then we have the regularity 1
(w22,Ωf + ζ22,Ωp ) 2 + r1,Ωf ≤ c(g0 + θ1,Ωf ).
(4.15)
Moreover, there are the following approximation properties: there exists (wh , ζh , rh ) ∈ Hf,h × Hp,h × Qh such that w − wh 1,Ωf ≤ chw2,Ωf ,
(4.16)
ζ − ζh 1,Ωp ≤ chζ2,Ωp ,
(4.17)
r − rh 0,Ωf ≤ chr1,Ωf .
(4.18)
Lemma 4.4. Under the assumption of regularity (4.15), we have the following estimate 1
(u − uh 20,Ωf + ϕ − ϕh 20,Ωp ) 2 + p − ph −1,Ωf ≤ ch2 .
(4.19)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1138
Proof. Comparing the problems (2.12)–(2.13) and (3.3)–(3.4), then taking (vh , qh ) = (wh , rh ) = (wh , ζh , rh ) gives a(e, wh ) + b(wh , p − ph ) = 0,
(4.20)
b(u − uh , rh ) = 0.
(4.21)
Setting v = e, q = p − ph in (4.14), we have (g, e) + (θ, p − ph )Ωf = a(e, w) − b(u − uh , r) + b(w, p − ph ) = a(e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph ).
(4.22)
The right side in (4.22) are bounded by Schwarz and trace inequalities, |a(e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph )| = |af (u − uh , w − wh ) + ap (ϕ − ϕh , ζ − ζh ) + aΓ (e, w − wh ) − b(u − uh , r − rh ) + b(w − wh , p − ph )| ≤ c(u − uh 1,Ωf (w − wh 1,Ωf + ζ − ζh 1,Ωp + r − rh 0,Ωf ) + ϕ − ϕh 1,Ωp (w − wh 1,Ωf + ζ − ζh 1,Ωp ) + p − ph 0,Ωf w − wh 1,Ωf ).
(4.23)
Combining (4.13), (4.16)–(4.18), (4.23) with (4.22) yields (g, e) + (θ, p − ph )Ωf ≤ ch2 (w2,Ωf + ζ2,Ωp + r1,Ωf ) 1
≤ ch2 ((w22,Ωf + ζ22,Ωp ) 2 + r1,Ωf ) ≤ ch2 (g0 + θ1,Ωf ),
(4.24)
(4.19) thus follows. 2 Next, we give the convergence of the modified two-grid scheme, it can be found in [39,38]. Lemma 4.5. Let (uh , ph , ϕh ) be the solution of modified two-grid scheme, then the following error estimates hold: u − uh 1,Ωf ≤ c(h + H 2 ),
(4.25)
p − ph 0,Ωf ≤ c(h + H 2 ),
(4.26)
ϕ − ϕh 1,Ωp ≤ c(h + H 2 ).
(4.27)
Finally, we will consider the convergence of the modified local and parallel scheme. Theorem 4.1. Let (uh , ph , ϕh ) be the solution of the modified local and parallel scheme, we have the following error estimates |u − uh |1,Ωf + |p − ph |0,Ωf ≤ c(h + H 2 ).
(4.28)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1139
Proof. To prove (4.28), we first prove: u − uh 1,Dj + p − ph 0,Dj ≤ c(h + H 2 ).
(4.29)
Let vh = (vh , 0) in (3.3), we have af (uh , vh ) + b(vh , ph ) = (f1 , vh )Ωf − g
ϕh vh · nf ,
∀vh ∈ Hf,h (Ωf )
(4.30)
∀qh ∈ Qh (Ωf ).
(4.31)
Γ
b(uh , qh ) = 0, Comparing the Algorithm 3 and (4.30)–(4.31), gives af (uh − uh , vh ) + b(vh , ph − ph ) = g
(ϕh − ϕh )vh · nf
∀vh ∈ Hf,h (Ωj ),
(4.32)
∀qh ∈ Qh (Ωj ).
(4.33)
ΓΩj
b(uh − uh , qh ) = 0 Note that from [30,39,38], we have ϕh − ϕh 1,Ωf ≤ cH 2 .
(4.34)
By Lemma 4.2, trace inequality (2.15), (4.34) and (4.19), we have uh − uh 1,Dj + ph − ph 0,Dj ≤ c(ϕh − ϕh L2 (ΓΩj ) + uh − uh 0,Ωj + ph − ph −1,Ωj ) ≤ c(ϕh − ϕh 1,Ωf + uh − uH 0,Ωj + ph − pH −1,Ωj + ejh 0,Ωj + ηhj −1,Ωj ) ≤ c(H 2 + ejh 0,Ωj + ηhj −1,Ωj ).
(4.35)
To estimate ejh 0,Ωj + ηhj −1,Ωj , we consider the following problem ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ∂w ⎪ −ντi ∂n ⎪ f ⎪ ⎪ ⎩
−νΔw − ∇r = θ in Ωj , ∇ · w = φ in Ωj , = g(ϕh − ϕH )w/v on ΓΩj ,
∂w r − νnf ∂n f
νg = α tr(K) w · τi 1 ≤ i ≤ (d − 1) on ΓΩj ,
(4.36)
w = 0 on ∂Ωj \ ΓΩj .
There holds w2,Ωj + r1,Ωj ≤ c(θ0,Ωj + φ1,Ωj ).
(4.37)
Then the weak form of (4.36) can be read as follows: for ∀(v, q) ∈ Hf (Ωj ) × L2 (Ωj ) af (w, v) − b(v, r) + b(w, q) + g ΓΩj
(ϕh − ϕH )w · nf = (θ, v)Ωj + (φ, q)Ωj .
(4.38)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1140
Let (wμ , rμ ) ∈ Hf,μ (Ωj ) × Qμ (Ωj ) be the finite element approximation of (w, r), i.e., af (w − wμ , v) − b(v, r − rμ ) + g
(ϕh − ϕH )(w − wμ ) · nf = 0 ∀v ∈ Hf,μ (Ωj ),
(4.39)
ΓΩj
b(w − wμ , q) = 0 ∀q ∈ Qμ (Ωj ),
(4.40)
w − wμ 1,Ωj + r − rμ 0,Ωj ≤ cμ(w2,Ωj + r1,Ωj ) ≤ cμ(θ0,Ωj + φ1,Ωj ).
(4.41)
where μ = h or H. Then
It follows from the problem (4.30)–(4.31) and (3.10)–(3.11) that af (ejh , v) + b(v, ηhj ) − b(ejh , q) = af (uh − uH , v) + b(v, ph − pH ) − b(uh − uH , q) +g (ϕh − ϕh )v · nf , ∀(v, q) ∈ Hf,h (Ωj ) × Qh (Ωj ), ΓΩj
(4.42)
af (uh − uH , wH ) + b(wH , ph − pH ) − b(uh − uH , rH ) + g
(ϕh − ϕH )wH · nf = 0.
(4.43)
ΓΩj
Setting v = ejh , q = ηhj in (4.38), then combining (4.39)–(4.40) with (4.42)–(4.43) yields (θ, ejh ) + (φ, ηhj ) =
af (ejh , w)
−
b(ejh , r)
+
b(w, ηhj )
(ϕh − ϕH )w · nf
+g ΓΩj
=
af (ejh , wh )
−
b(ejh , rh )
+
b(wh , ηhj )
+g
(ϕh − ϕH )wh · nf
ΓΩj
= af (uh − uH , wh ) + b(wh , ph − pH ) − b(uh − uH , rh ) + g
(2ϕh − ϕh − ϕH )wh · nf
ΓΩj
= af (uh − uH , wh − wH ) + b(wh − wH , ph − pH ) − b(uh − uH , rh − rH ) +g (ϕh − ϕH )(wh − wH ) · nf + g (ϕh − ϕh )wh · nf . ΓΩj
(4.44)
ΓΩj
Hence, we can derive from (4.41), (4.44), trace inequality (2.15), (4.13) and (4.34) that |(θ, ejh ) + (φ, ηhj )| ≤ c[H(uh − uH 1,Ωf + ph − pH 0,Ωf + ϕh − ϕH 1,Ωf ) + ϕh − ϕh 1,Ωf ](θ0,Ωj + φ1,Ωj ) ≤ cH 2 (θ0,Ωj + φ1,Ωj ),
(4.45)
which implies ejh 0,Ωj + ηhj −1,Ωj ≤ cH 2 .
(4.46)
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1141
Applying (4.46) to (4.35), we get uh − uh 1,Dj + ph − ph 0,Dj ≤ cH 2 .
(4.47)
Then (4.29) follows from triangle inequalities and (4.13). Finally, from the definition of |u − uh |1,Ωf and |p − ph |0,Ωf in (4.12), there holds |u − uh |1,Ωf + |p − ph |0,Ωf ≤
m (u − uh 1,Dj + p − ph 0,Dj ) j=1
≤ c(h + H 2 ).
2
(4.48)
5. Numerical experiments In this section, we give three numerical experiments. In the first two experiments, we consider the cases for which the exact solution is known and unknown respectively, to demonstrate the accuracy and efficiency of the modified local and parallel method. In the final case, we choose different K to examine the effects of the hydraulic conductivity tensor. Assume the computational domain Ωf = (0, 1) × (1, 2) and Ωp = (0, 1) × (0, 1) with interface Γ = (0, 1)×{1}. To apply the modified local and parallel algorithm, we divide the domain Ωf into four subdomains D1 = (0, 1/4) × (1, 2),
D2 = (1/4, 1/2) × (1, 2),
D3 = (1/2, 3/4) × (1, 2),
D4 = (3/4, 1) × (1, 2),
and set Ω1 = (0, 5/16) × (1, 2),
Ω2 = (3/16, 9/16) × (1, 2),
Ω3 = (7/16, 13/16) × (1, 2),
Ω4 = (11/16, 1) × (1, 2).
The well-known Mini elements (P 1b − P 1) are used for the Stokes problem and the linear Lagrangian elements (P 1) for the Darcy flow. The relation of the coarse mesh H and the fine mesh h is chosen as h = H 2 . Moreover, we apply the software package FreeFEM++ to implement the algorithms and all results are generated by the same computer. For simplicity of notation, we denote ehϑ = ϑh − ϑ,
h
ehϑ = ϑ − ϑ,
where ϑ can be ϕ, u, p. Experiment 1. We set the physical parameters g, ν, α equal to 1, and K = I. The exact solution is given by ⎧ 2 2 2 3 ⎪ ⎨ u = ([x (y − 1) + y], [− 3 x(y − 1) + 2 − π sin(πx)]), p = [2 − π sin(πx)] sin( π2 y), ⎪ ⎩ ϕ = [2 − π sin(πx)][1 − y − cos(πy)]. The forcing terms and the boundary conditions follow the exact solution.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1142
Table 1 The convergence performance and CPU time of the coupling scheme. h
∇(ϕ − ϕh )0,Ωp
Rate
∇(u − uh )0,Ωf
Rate
p − ph 0,Ωf
Rate
CPU
1 16 1 36 1 64
0.375149 0.163785 0.092593
– 1.02 0.99
0.360524 0.162884 0.091482
– 0.98 1.00
0.148571 0.055708 0.033513
– 1.21 0.88
5.893 34.106 243.437
Table 2 The convergence performance and CPU time of the modified two-grid method. h
h ∇eϕ 0,Ωp
Rate
∇eh u 0,Ωf
Rate
eh p 0,Ωf
Rate
Stokes
Darcy
CPU
1 16 1 36 1 64 1 100
0.375473 0.166299 0.094427 0.060576
– 1.00 0.98 1.02
0.370676 0.161416 0.091266 0.058505
– 1.03 0.99 1.02
0.166929 0.056210 0.032502 0.019653
– 1.34 0.95 1.15
0.125 0.655 2.262 5.975
0.015 0.078 0.265 0.655
0.281 1.279 3.978 9.922
Table 3 The convergence performance and CPU time of the modified local and parallel scheme. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
Stokes
Darcy
CPU
1 16 1 36 1 64 1 100
0.375473 0.166299 0.094427 0.060576
− 1.00 0.98 1.02
0.370912 0.163275 0.091858 0.059705
− 1.01 1.00 0.99
0.159621 0.060694 0.029189 0.017481
− 1.19 1.27 1.18
0.062 0.25 0.78 2.044
0.015 0.093 0.266 0.655
0.232 0.89 2.511 6.038
In Table 1, we list the errors between the exact solution and the finite element approximate solutions of the coupled problem, convergent rates and computational time with the mesh size h = 1/16, 1/36, 1/64. It’s obvious that the errors for ∇ϕ, ∇u and p in the L2 -norm are all around of the order of O(h). In Table 2, the errors and convergence rates of the modified two-grid method are given. Observe that, the numerical results exhibit the first-order reduction. Meanwhile, the computational time of the Stokes problem (3.8)–(3.9), the Darcy problem (3.7), and the Stokes–Darcy problem (3.5)–(3.9) are listed. It is showed that the Stokes problem takes more time than the Darcy problem. In Table 3, we show the convergence performance and computational time of the modified local and parallel scheme with the mesh size h = 1/16, 1/36, 1/64, 1/100. Especially, “Stokes” denotes the max CPU time spent by each local Stokes sub-problem (3.10)–(3.11), and “Darcy” means the computational time of the problem (3.7). Apparently, the numerical results agree with the theoretical analysis strongly. That is to say, the modified local and parallel scheme has the same accuracy as the coupling scheme. From Table 2 and Table 3, it’s obvious that the modified local and parallel scheme spends less time than the modified two-grid scheme. Experiment 2. In this numerical experiment, we consider a case for which the exact solution for the coupled problem is unknown. The physical parameters are the same as in Experiment 1, the boundary conditions are homogeneous Dirichlet boundary conditions: u=0
on Γf ,
ϕ = 0 on Γp .
The external forces are f1 = (−(A1 + B1 ) + C1 , −(A2 + B2 ) + C2 ), f2 = −12(2x − 1)y, where
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1143
Table 4 The convergence performance of the modified local and parallel scheme. h
∇(ϕh0 − ϕh )0,Ωp
Rate
|∇(uh0 − uh )|0,Ωf
Rate
|ph0 − ph |0,Ωf
Rate
1 16 1 25 1 36 1 49 1 64
0.076086 0.047841 0.035011 0.029221 0.020406
– 1.04 1.01 0.96 1.00
0.070361 0.046777 0.034825 0.027308 0.020989
– 0.95 0.92 0.91 0.92
0.022657 0.017673 0.012593 0.010220 0.007295
– 0.99 0.98 0.94 0.97
Table 5 The convergence performance of the modified local and parallel algorithm with k = 2. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
0.017839 0.007818 0.004467 0.002838
– 1.02 0.97 1.02
0.045087 0.020268 0.011283 0.007215
– 0.99 1.02 1.00
0.016630 0.006532 0.003287 0.001909
– 1.15 1.19 1.22
Table 6 The convergence performance of the modified local and parallel algorithm with k = 0.2. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
0.178335 0.078160 0.044599 0.028389
– 1.02 0.98 1.01
0.045089 0.020269 0.011284 0.007215
– 0.99 1.02 1.00
0.021258 0.009331 0.004455 0.002741
– 1.02 1.28 1.09
A1 = 2(6x2 − 6x + 1)(2y 3 − 6y 2 + 5y − 2), B1 = 12x2 (x − 1)2 (y − 1), C1 = 2(6x2 − 6x + 1)y, A2 = −6(2x − 1)(y − 2)2 (y 2 + 1), B2 = −2x(x − 1)(2x − 1)(6y 2 − 12y + 5), C2 = 2x(x − 1)(2x − 1). In this case, we first solve the coupled Stokes–Darcy model in a fine mesh grid h0 (h0 = 1/64), derive the approximate solution (uh0 , ph0 , ϕh0 ), and regard the approximate solution as one “exact solution”. In Table 4, we show the errors between the “exact solution” and the modified local and parallel scheme solutions and give the convergence rate. Experiment 3. The final case is carried out in order to assess the influence of the physical parameters on the convergence rate of the modified local and parallel scheme. In this case, we take K = kI and choose the exact solution ⎧ 2 ⎪ ⎨ u = ((y − 1) , x(x − 1)), p = 2ν(x + y − 1), ⎪ ⎩ ϕ = k−1 (x(1 − x)(y − 1) + (y − 1)3 /3) + 2ν x. g The forcing terms and the boundary conditions follow the exact solution, the physical parameters g, ν, α are the same as in Experiment 1. In Tables 5–7, we list the convergence performance of the modified local and parallel algorithm with K = 2I, K = 0.2I, K = 0.02I, respectively. The errors of ∇ϕ and p in L2 -normal increase as k decreases,
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1144
Table 7 The convergence performance of the modified local and parallel algorithm with k = 0.02. h
∇eh ϕ 0,Ωp
Rate
|∇eh u |0,Ωf
Rate
|eh p |0,Ωf
Rate
1 16 1 36 1 64 1 100
1.783670 0.781297 0.446450 0.284395
– 1.02 0.97 1.01
0.045188 0.020327 0.011315 0.007235
– 0.99 1.02 1.00
0.112322 0.056715 0.030252 0.019952
– 0.84 1.09 0.93
while the errors of ∇u in L2 -normal almost remain the same. The changes of k do not influence the convergence rate of ∇ϕ, p and ∇u. On the whole, we proposed a modified local and parallel algorithm for the mixed Stokes and Darcy problem. Optimal error estimates are derived. Numerical experiments show that the modified local and parallel scheme we proposed not only has the same accuracy as the coupling scheme, but also saves a lot of time than the modified two-grid scheme. We thus conclude that the modified local and parallel algorithm is very effective. References [1] T. Arbogast, D.S. Brunson, A computational method for approximating a Darcy–Stokes system governing a vuggy porous medium, Comput. Geosci. 11 (2007) 207–218. [2] G. Beavers, D. Joseph, Boundary conditions at a naturally impermeable, J. Fluid Mech. 30 (1967) 197–207. [3] Y. Boubendir, S. Tlupova, Domain decomposition methods for solving Stokes–Darcy problems with boundary integrals, SIAM J. Sci. Comput. 35 (2013) B82–B106. [4] M.C. Cai, M. Mu, A multilevel decoupled method for a mixed Stokes/Darcy model, J. Comput. Appl. Math. 236 (2012) 2452–2465. [5] M.C. Cai, M. Mu, J.C. Xu, Preconditioning techniques for a mixed Stokes/Darcy model in porous media applications, J. Comput. Appl. Math. 233 (2009) 346–355. [6] M.C. Cai, M. Mu, J.C. Xu, Numerical solution to a mixed Navier–Stokes/Darcy model by the two-grid approach, SIAM J. Numer. Anal. 47 (2009) 3325–3338. [7] Y. Cao, M. Gunzburger, X. Hu, F. Hua, X. Wang, W. Zhao, Finite element approximation for Stokes–Darcy flow with Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 47 (2010) 4239–4256. [8] Y. Cao, M. Gunzburger, F. Hua, X. Wang, Coupled Stokes–Darcy model with Beavers–Joseph interface boundary condition, Commun. Math. Sci. 8 (2010) 1–25. [9] W. Chen, M. Gunzburger, F. Hua, X.M. Wang, A parallel Robin–Robin domain decomposition method for the Stokes– Darcy system, SIAM J. Numer. Anal. 49 (2011) 1064–1084. [10] P. Chidyagwai, B. Rivière, A two-grid method for coupled free flow with porous media flow, Adv. Water Resour. 34 (2011) 1113–1123. [11] M. Discaaaiati, A. Quarteroni, A. Valli, Robin–Robin domain decomposition methods for the Stokes–Darcy coupling, SIAM J. Numer. Anal. 45 (2007) 1246–1268. [12] M. Discacciati, Domain decomposition methods for the coupling of surface and groundwater flows, Ph.D. dissertation, École Polytechnique Fédérale de Lausanne, 2004. [13] M. Discacciati, A. Quarteroni, Convergence analysis of a subdomain iterative method for the finite element approximation of the coupling of Stokes and Darcy equations, Comput. Vis. Sci. 6 (2004) 93–103. [14] M. Discacciati, E. Miglio, A. Quarteroni, Mathematical and numerical models for coupling surface and groundwater flows, Appl. Numer. Math. 43 (2002) 57–74. [15] M. Discacciati, A. Quarteroni, Analysis of a domain decomposition method for the coupling Stokes and Darcy equations, in: F. Brezzi, et al. (Eds.), Numerical Analysis and Advanced Applications, Enumath 2001, Springer, Milan, 2003, pp. 3–20. [16] J. Galvis, M. Sarkis, Balancing domain decomposition methods for mortar coupling Stokes–Darcy Systems, in: Proceedings of the 16th International Conference on Domain Decomposition Methods, New York, 2005. [17] V. Girault, B. Rivière, DG approximation of coupled Navier–Stokes and Darcy equations by Beaver–Joseph–Saffman interface condition, SIAM J. Numer. Anal. 47 (2009) 2052–2089. [18] R. Glowinski, T. Pan, J. Periaux, A Lagrange multiplier/fictitious domain method for the numerical simulation of incompressible viscous flow around moving grid bodies: I. Case where the rigid body motions are known a priori, C. R. Acad. Sci. Paris Ser. I Math. 324 (1997) 361–369. [19] N. Hanspal, A. Waghode, V. Nassehi, R. Wakeman, Numerical analysis of coupled Stokes/Darcy flow in industrial filtrations, Transp. Porous Media 64 (2006) 73–101. [20] Y.N. He, J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms for the Navier–Stokes problem, J. Comput. Math. 24 (2006) 227–238. [21] Y.N. He, J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms for the Stokes problem, Numer. Math. 109 (2008) 415–434. [22] F. Hua, Modeling, analysis and simulation of Stokes–Darcy system with Beavers–Joseph interface condition, Ph.D. dissertation, The Florida State University, 2009.
G. Du et al. / J. Math. Anal. Appl. 435 (2016) 1129–1145
1145
[23] B. Jiang, A parallel domain decomposition method for coupling of surface and groundwarter flows, Comput. Methods Appl. Mech. Engrg. 198 (2009) 947–957. [24] W.J. Layton, F. Schieweck, I. Yotov, Coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 40 (2003) 2195–2218. [25] Q.F. Liu, Y.R. Hou, Local and parallel finite element algorithms for time-dependent convection–diffusion equations, Appl. Math. Mech. (English Ed.) 30 (2009) 787–794. [26] F. Ma, Y. Ma, W. Wo, Local and parallel finite element algorithms based on two-grid discretization for steady Navier–Stokes equations, Appl. Math. Mech. (English Ed.) 28 (2007) 27–35. [27] S. Markus, E. Houstis, A. Catlin, J. Rice, P. Tsompanopoulou, E. Vavalis, D. Gottfried, K. Su, G. Balakrishnan, An agent-based netcentric framework for multidisciplinary problem solving environments (MPSE), Int. J. Comput. Eng. Sci. 1 (2000) 33–60. [28] A. Marquez, S. Meddahi, F.J. Sayas, A decoupled preconditioning technique for a mixed Stokes–Darcy model, J. Sci. Comput. 57 (2013) 174–192. [29] M. Mu, Solving composite problems with interface relaxation, SIAM J. Sci. Comput. 20 (1999) 1394–1416. [30] M. Mu, J.C. Xu, A two-grid method of a mixed Stokes–Darcy model for coupling fluid flow with porous media flow, SIAM J. Numer. Anal. 45 (2007) 1801–1813. [31] M. Mu, X.H. Zhu, Decoupled schemes for a non-stationary mixed Stokes–Darcy model, Math. Comp. 79 (2010) 707–731. [32] P. Saffman, On the boundary condition at the surface of a porous media, Stud. Appl. Math. 50 (1971) 93–101. [33] L. Shan, H.B. Zheng, Partitioned time stepping method for fully evolutionary Stokes–Darcy flow with the Beavers–Joseph interface conditions, SIAM J. Numer. Anal. 51 (2013) 813–839. [34] L. Shan, H.B. Zheng, W.J. Layton, A decoupling method with different subdomain time steps for the nonstationary Stokes–Darcy model, Numer. Methods Partial Differential Equations 29 (2013) 549–583. [35] J.C. Xu, A.H. Zhou, Some local and parallel properties of finite element discretizations, in: Proceedings for Eleventh International Conference on Domain Decomposition Methods, Greenwich, 1999, pp. 140–147. [36] J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms based on two-grid discretizations, Math. Comp. 69 (1999) 881–909. [37] J.C. Xu, A.H. Zhou, Local and parallel finite element algorithms based on two-grid discretizations for nonlinear problems, Adv. Comput. Math. 14 (2001) 293–327. [38] T. Zhang, J.Y. Yuan, Two novel decoupling algorithms for the steady Stokes–Darcy model based on two-grid discretizations, Discrete Contin. Dyn. Syst. Ser. B 19 (2014) 849–865. [39] L.Y. Zuo, Y.R. Hou, A decoupling two-grid algorithm for the mixed Stokes–Darcy model with the Beavers–Joseph interface condition, Numer. Methods Partial Differential Equations 3 (2014) 1066–1082.