- Email: [email protected]
New asymptotic expansion for the diffusion problem with high-contrast inclusions: Numerical results
Journal of Computational and Applied Mathematics 365 (2020) 112363
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
New asymptotic expansion for the diffusion problem with high-contrast inclusions: Numerical results ∗
V.K. Kramarenko a , , Yu. A. Kuznetsov b a b
Marchuk Institute of Numerical Mathematics, Russian Academy of Science, Gubkina 8, Moscow, Russian Federation University of Houston, Department of Mathematics, Houston, USA
article
info
Article history: Received 27 February 2019 Received in revised form 9 July 2019 Keywords: Diffusion equation High-contrast inclusions Finite elements Asymptotic expansions Saddle-point matrix
a b s t r a c t This paper is concerned with the new approach to the numerical solution of diffusion equations with high contrast coefficients. The classical P1 finite element system of equations is replaced by the equivalent system of the saddle point type. For the saddle point problem a new recently invented asymptotic expansion with respect to a small parameter is discussed and the estimates for the convergence factor are given. Numerical results completely confirm the theoretical accuracy and convergence results for the considered method. © 2019 Elsevier B.V. All rights reserved.
0. Introduction Elliptic problems with high contrast coefficients recently have been a topic of active research, see [1] and references therein. An approach based on the classical asymptotic expansion technique is among the most popular methods for solving this type of problems. A new approach to the asymptotic expansion for the P1 finite element approximations of the diffusion problem with high contrast inclusions was proposed and investigated in [2]. In what follows the main idea is to replace the classical P1 finite element problem by the equivalent problem of the saddle point type. This idea was invented originally for solving of algebraic systems with singularly perturbed symmetric positive definite matrices in [3]. In [4], this procedure was applied to solving P1 finite element discretizations of the diffusion equations with high contrast inclusions. An efficient preconditioned Lanczos method was proposed and investigated. Rigorous analysis with new convergence proof of the preconditioned iterative methods for underlying algebraic saddle point systems of equations was appeared in [5]. The current paper complements paper [2]. By virtue of the fact it contains the description of some practical aspects invented in [2] and numerical results which completely confirm the earlier theoretical results. The paper is organized as follows. In Section 1 we formulate the problem and discuss the algebraic aspects of the replacement of the original P1 finite element problem by the equivalent saddle type problem. In Section 2 we recall the new asymptotic expansion proposed in [2] and the easier convergence results. We discuss the estimates of the convergence factor for the new asymptotic expansion in Section 3. The model test problem and the numerical results are given in Section 4. For the sake of simplicity we consider the case of a square domain with multiple uniformly and randomly distributed square inclusions. The numerical results confirm the theoretical statements and the proposed estimates for the convergence factor. ∗ Corresponding author. E-mail addresses: [email protected] (V.K. Kramarenko), [email protected] (Yu. A. Kuznetsov). https://doi.org/10.1016/j.cam.2019.112363 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
1. Problem formulation We consider the diffusion problem
−∇ [k (x) ∇ u] = f , x ∈ Ω , u = 0, x ∈ ∂Ω, where Ω ∈ R2 is a polygonal domain with the boundary ∂ Ω and k (x) ⩾ 1 is a positive piecewise constant function. We introduce polygonal subdomains (inclusions) ωs of Ω , s = 1, m, where m is a positive integer, such that ωs ∩ωt = ∅ and ωs ∩ ∂ Ω = ∅, s, t = 1, m. For the sake of simplicity we assume that subdomains ωs are convex and
{ k (x) =
ks = 1 +
1,
1
εs
≡ const > 1 in ωs , s = 1, m,
otherwise,
and consider the weak formulation of : Find u ∈ H01 (Ω ) such that
∫ Ω
(k (x) ∇ u) · ∇v dx ≡
∫ ∇ u · ∇v dx + Ω
∫ ∫ m ∑ 1 ∇ u · ∇v dx = f v dx εs ωs Ω
(1)
s=1
∀u, v ∈ H01 (Ω ). Let Ωh be a triangular mesh in Ω conforming with the boundary ∂ωs , i.e. ∂ωs is the union of edges from Ωh , s = 1, m. We denote by Vh the classical P1-finite element subspace of H01 (Ω ) on Ωh . The P1-finite element method: Find uh ∈ Vh , such that
∫ Ω
(k (x) ∇ u) · ∇v =
∫
f v dx
Ω
∀v ∈ Vh
results in the algebraic system Aε u = f
(2)
with the N × N matrix Aε = Aε > 0 and the vector f ∈ R , where N is the number of interior mesh nodes in Ωh . With an appropriate ordering the matrix Aε can be presented as 2 × 2 block matrix T
[ Aε =
A11,ε A21
A12 A22
N
]
with n × n submatrix A11,ε = A11 + B1,ε , where B1,ε is the m × m block diagonal matrix with the diagonal blocks 1ε As ∈ Rns ×ns , s = 1, m. Here, ns is the number of ∑m nodes belonging to ωs , s = 1, m, n = s=1 ns , and matrices As are defined by the identities
(As us , v s ) =
∫ ωs
∇ uh · ∇vh dx ∀us , v s ∈ Rns ,
(3)
where uh , vh ∈ Vs,h and Vs,h is the restriction of Vh onto ωs , s = 1, m. In other words,
[ Aε =
]
A11 A21
[
A12 B + 1,ε A22 0
]
0 , 0
(4)
where
{ B1,ε = diag
1
ε1
A1 , . . . ,
1
εm
}
Am .
Using the procedure proposed in [3,4], we replace (2) by the system Au Bu
+ −
BT p Σε B1 p
= =
f, 0,
(5)
where
[
A11 A= A21
A12 A22
] (6)
is n × n matrix, B = [B1 0]
(7)
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
3
is n × N matrix, B1 = diag {A1 , . . . , Am } ,
Σε = diag {ε1 I1 , . . . , εm Im } , and p = Σε−1 0 u.
]
[
Here, Is are the identity ns × ns matrices, s = 1, m. It is obvious, that the matrix
(
A B
A=
)
BT
−Σε B1
is singular, dim(ker A) = m, ker BT = ker B1 , and any vector w ∈ ker B1 can be presented in the block form by
⎤ w1 ⎢ ⎥ w = ⎣ ... ⎦ , wm ⎡
where w s ∈ ker As are vectors with constant components, s = 1, m. The solution vector u ∈ Rn in (5) is unique, and the solution vector p ∈ Rn in (5) is unique up to an arbitrary additive vector w ∈ ker B1 . Henceforth, we shall always assume that the solution vector p in (5) is orthogonal to ker B1 , i.e. the uniqueness of p. Under the latter assumption p ⊥ ker B1 , system (5) is equivalent to the system
[ ] Lα
u ≡ p
(
BT −Σε B1 − Qα
A B
)[ ]
[ ]
u f , = p 0
(8)
where
{
}
Qα = diag α1 Q1 , . . . , αm Qm ,
αs = consts > 0, s = 1, m, Qs = Ms w s w Ts Ms ,
s = 1 , m,
and Ms are the underlying mass matrices defined by
(Ms us , v s ) =
∫ ωs
us,h vs,h dx
∀us,h , vs,h ∈ Vs,h ,
s = 1, m. We provide the following alternative definition of system (8). We replace the variational problem (1) by the equivalent variational saddle point problem: find u ∈ H01 (Ω ), ps ∈ H 1 (ωs ), s = 1, m, such that
∫ ∇ u · ∇v dx + Ω
∫ ωs
m ∫ ∑
ωs
s=1
∇ u · ∇ qs dx − εs
∫
∫ ωs
∇ ps · ∇v dx =
Ω
∇ ps · ∇ qs dx − αs
f v dx,
(9)
∫
∫ ωs
ps dx ·
ωs
qs dx = 0
(10)
∀v ∈ H01 (Ω ), qs ∈ H 1 (ωs ), s = 1, m, where αs are positive constants, s = 1, m.
Then system (8) results from the finite element discretization of (9)–(10): find uh ∈ Vh , ps,h ∈ Vs,h , s = 1, m, such that
∫ Ω
∫ ωs
∇ uh · ∇vh dx +
m ∫ ∑
ωs
s=1
∇ uh · ∇ qs,h dx − εs
∫
∫ ωs
∇ ps,h · ∇vh dx =
Ω
f vh dx,
∇ ps,h · ∇ qs,h dx − αs
∫
∫ ωs
ps,h dx ·
ωs
qs,h dx = 0
∀vh ∈ Vh , qs,h ∈ Vs,h , s = 1, m. 2. New asymptotic expansion 1 The matrix Lα in (8) is nonsingular for any εs ⩾ 0, and the entries of L− α are infinitely differentiable functions of εs ⩾ 0 1 and tend to the entries of L− as ε → 0, s = 1 , m. To this end, it is natural to use the classical asymptotic expansion for 0
4
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
the solution vectors u and p in (8) (the solution finite element functions uh and ps,h , s = 1, m, in (9)–(10)) with respect to parameters εs , s = 1, m. In this section, similar to [2], we consider only the case of ε1 = ε2 = · · · = εm = ε . First, we consider the system (8) with ε1 = ε2 = · · · = εm = 0, i.e. the system 0
Au 0 Bu
0
BT p 0 Qα p
+ −
f, 0,
= =
(11)
0
0
0
and estimate the errors u − u and p − p in A1 -norm and B1 -seminorm, respectively. We recall that Qα p = 0. Below we present some results and proofs from [2]. The error vectors satisfy the system
(
)(
BT
A B
−Σε B1
0
u−u 0 p−p
)
( =
0 Σε B1 p0
)
.
(12)
0
Eliminating the vector u − u in (12) we obtain the system
= −Σε B1 p0 ,
(S + Σε B1 ) p − p
) 0
(
(13)
where −1 S = BA−1 BT ≡ B1 S11 B1
and 1 S11 = A11 − A12 A− 22 A21 .
Consider the eigenvalue problem −1 B1 S11 B1 w = µB1 w
(14)
with the eigenvalues 0 = µ1 = µ2 = · · · = µm < µm+1 ⩽ 1 where µm+1 is the minimal positive eigenvalue. The inequality µm+1 ⩽ 1 was proved in [3]. Then from (13) we derive the estimate
p − p0
B1
⩽
ε
· p0 B ,
µm+1
(15)
1
where ∥·∥B1 denotes 0 the semi-norm generated by the matrix B1 . 0 To estimate p B in (15), we return to system (11). Eliminating the vector u we have the equation 1
0
−1
Sp = −B1 A
f.
Then the estimate
0 p
B1
=
1
µm+1
f
(16)
A− 1
follows from the fact that
1 12 − 1 B A 2 = A− 21 B 2 ⩽ 1 1 1 2
(17)
2
where ∥·∥A−1 is the norm generated by A−1 . Thus, we obtain the estimate
p − p0
B1
⩽
ε · f A−1 . µ2m+1 0
Multiplying the first equation in (12) by the vector u − u , and applying the Schwarz inequality and (17), we get the estimates
u − u0 ⩽ p − p 0 ⩽ A B 1
ε µ
2 m+1
· f A− 1 .
(18)
We now consider the expansions u=
∞ ∑
ε i u(i) ,
(19)
ε i p(i) ,
(20)
i=0
p=
∞ ∑ i=0
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 (i)
where the vectors u , p (i)
Au (i) Bu
+
T
B p
(i)
(i)
5
satisfy Eq. (11) for i = 0 and the equations
= =
0, (i−1) B1 p ,
i = 1, 2, . . .. Similar to above we readily derive the estimates
(i) p
B1
⩽
1
µm+1
· p(i−1) B ⩽ 1
1
(µm+1 )i
· p(0) B , 1
(21)
i = 1, 2, . . .. Define the error vectors (r)
y
=u−
r ∑
ε l u(l) ,
l=0
z
(r)
=p−
r ∑
ε l p(l)
l=0
for some integer r ⩾ 1. (r) (r) It is obvious that the vectors y and z satisfy the system (r)
Ay (r) By
+ −
(r)
BT z ε B1 z (r)
= =
0, −εr +1 B1 p(r) ,
(22)
r = 1, 2, . . .. (r) Again eliminating the vector y in (22) we obtain the system
(S + ε B1 ) z (r) = −ε r +1 B1 p(r) . Now, using the estimate
(r) z
B1
⩽
ε r +1 · p(r) B , 1 µm+1
and (16) and (21), and the estimate
(r) y ⩽ z (r) , A B
(23)
1
we obtain the final estimates r +1 r y ⩽ z (r) ⩽ ε · f A−1 , r +2 A B1 µm+1
r = 0, 1, . . . .
(24)
On the bases of estimates (18) and (24) we reach the following conclusion: Statement 1. The estimate r +1 (r) y ⩽ ε · u∗ A r +2 A µm+1
(25)
∗
holds for any r ⩾ 0 where u = A−1 f . Here, u∗h is the P1-finite element solution of the Poisson problem
− △ u∗ = f in Ω , u∗ = 0 on ∂ Ω on the mesh Ωh . Statement 2. The condition
ε
µm+1
<1
(26) (r)
is sufficient for convergence of uh to uh as r −→ +∞. 3. Estimation from below for µm+1 We replace (14) by an equivalent eigenvalue problem: Find µ ∈ R, v ∈ RN , w ∈ Rn , such that Av Bv
+
BT w
= =
0, −µB1 w,
(27)
6
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
where
v = −A−1 BT w ≡
[ ] v1 v2
−1 and v 1 = −S11 B1 w . It follows that
A21 v 1 + A22 v 2 = 0, i.e., v2,h ∈ VΩ \ω,h is the h-harmonic extension of v1,h ∈ Vω,h . Here Vω,h and VΩ \ω,h are the restrictions of Vh onto ω and Ω \ ω, respectively. The finite element formulation of (27) is as follows: Find µ ∈ R, vh ∈ Vh and wh ∈ Vω,h , such that
∫
∫
∫Ω ω
∇vh · ∇ϕh dx
+ ω
∇wh · ∇ϕh dx
0,
=
(28)
∫ ∇vh · ∇ψh dx
=
−µ ω
∇wh · ∇ψh dx
for any ϕh ∈ Vh , ψh ∈ Vω,h . To consider only nonzero eigenvalues, we impose additional conditions
∫
wh dx = 0,
ωs
s = 1, m.
(29)
It is clear that µm+1 in (14) is the minimal eigenvalue in (28)–(29). It was established in [4,5] with the reference to the finite element norm preserving extension result in [6] that in the case of regular shaped quasiuniform meshes Ωh the value of µm+1 is bounded form below by a positive constant independent of the mesh. Moreover, it was further established that under the additional assumptions that ωs are convex and regular shaped, and ds,t = inf |x − y| ⩾ αs ds ,
(30)
x∈ωs , y∈ωt , t ̸ =s
where αs are the positive constants and ds = |ωs |1/2 , s, t = 1, m, the value of µm+1 is bounded from below by a positive constant also independent of ds , s = 1, m. An alternative proof of the above results was recently proposed in [2]. In this section we follow the latest proof. Let µ be an eigenvalue and vh , wh corresponding eigenfunctions in (28)–(29). If we choose ϕh = vh and ψh = wh , we obtain the equality
∫
|∇vh |2 dx = µ
Ω
∫ ω
|∇wh |2 dx.
Taking ϕh = vh and using the inequality
⏐ [∫ ⏐∫ ]1/2 [∫ ]1/2 ⏐ ⏐ 2 2 ⏐ ∇wh · ∇vh dx⏐ ⩽ |∇w | |∇v | dx dx h h ⏐ ⏐ ω
ω
Ω
we obtain
∫ Ω
|∇vh |2 dx ⩽
∫ ω
|∇wh |2 dx
and conclude that the eigenvalues µ in (28)–(29) belong to the segment (0; 1]. Now we choose ϕh ∈ Vh such that ϕh = wh in ω and ϕh = w2,h in Ω \ ω, where w2,h is a finite element extension of wh from ω into Ω \ ω, i.e. w2,h = wh on ∂ω and w2,h ∈ VΩ \ω,h . From the equality
∫
∫ Ω
∇vh · ∇ϕh dx = −
ω
|∇wh |2 dx,
we derive the estimate
∫ ω
[ ] |∇wh |2 dx ⩽ 1 + C 2
∫ Ω
|∇vh |2 dx
with
∫ C2 =
Ω \ω
∫ ω
⏐ ⏐ ⏐∇w2,h ⏐2 dx . |∇wh |2 dx
(31)
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
7
Fig. 1. An example of ωs and ω ˆ s.
The value of C 2 is minimal when w2,h is the h-harmonic finite element extension of wh from ω to Ω \ ω. Choosing φh = vh and ψh = wh in (28) and using (31) we obtain the estimate 1
µm+1
⩽ 1 + C 2.
To estimate C 2 from above, we construct an appropriate function w2,h . We define in Ω a set of nonoverlapping open polygons ω ˆ s , s = 1, m, with the boundaries ∂ ωˆ s , such that ωs ∈ ωˆ s , i.e., ∂ωs ∩ ∂ ωˆ s = 0 and ∂ ωˆ s ∩ ∂ Ω = 0, s = 1, m. We assume that ω ˆ s are conforming with the mesh Ωh , i.e., ∂ ωˆ s are the unions of sides of triangular mesh cells in Ωh , s = 1, m. An example of hexagonal ωs and ω ˆ s , 1 ⩽ s ⩽ m, is given in Fig. 1, where the boundary ∂ ωˆ s is denoted by dash lines. We subject w2,h to the conditions
w2,h = 0 in Ω \ ωˆ ⋃ where ω ˆ = m ˆ s, s=1 ω w2,h = wh on ∂ωs , s = 1, m (s)
and w2,h ≡ w2,h in ω ˆ s , where
∫ ωˆ s
∇w2(s),h · ∇vh dx = 0 ∀vh ∈ Vωˆ s \ωs ,
(s) 2,h
ˆ s , s = 1, m. i.e., w are the h-harmonic extensions of wh from ωs into ω We additionally assume that ∂ ω ˆ s satisfy the conditions β · ds ⩽ inf |x − y| x∈∂ωs y∈∂ ω ˆs
with positive constants β independent of Ωh , and the meshes ω ˆ s \ ωs are regular shaped and quasiuniform, s = 1, m. Under these assumptions, it was shown in [2,4,5] that the values of 2
∫ cs2
=
ωˆ s \ωs
∫
ωs
|∇w2(s),h | dx |∇wh |2 dx
,
s = 1, m, are bounded from above by positive constants independent of the mesh Ωh and values of ds , s = 1, m. It follows immediately that 1
µm+1
⩽ 1 + max cs2 .
(32)
1⩽s⩽m
The values of 1/(1 + cs2 ) are the minimal eigenvalues of the eigenvalue problems: Find νs ∈ R, vs,h ∈ Vωˆ s ,h , ws,h ∈ Vωs ,h , such that
∫ ∫ ωˆ s ∇vs,h · ∇φh dx ∇vs,h · ∇ψh dx ωˆ s
+
∫
ωs
∇ws,h · ∇φh dx
= =
−νs
0, ∇w s,h · ∇ψs,h dx ω
∫
(33)
8
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
Fig. 2. An example of uniform distribution of ωs , s = 1, m, with d = 2H =
1 , 8
h=H=
1 , 16
m = 16.
∀φh ∈ Vωˆ s ,h , ψh ∈ Vωs ,h under the condition ∫ ws,h dx = 0,
(34)
ωs
s = 1, m. The eigenvalue problem (33), (34) is equivalent to the algebraic eigenvalue problem A(s) us B(s) v s
+
[
B(s)
]T
ws
= =
0, −νs As v s ,
(35)
under the condition (34), where As are defined in (3) and A(s) and B(s) are the corresponding blocks of the matrices A and B in (6) and (7), respectively, s = 1, m. For instance, B(s) = As
[
0 ,
]
s = 1, m. The minimal eigenvalues in (33)–(35) can be computed numerically, or estimated, s = 1, m. 4. Numerical results 4.1. Model test problem Let Ω be the unit square and ΩH be the square mesh with the mesh size H =
1 NH
, where NH is a positive integer. We
define the set of ωs , s = 1, m, assuming that each ωs is a square and a union of mesh cells in ΩH , 1 ⩽ s ⩽ m. In numerical tests we always assume that ωs have the same size d ≡ ds = kH, s = 1, m, where k ⩾ 2 is a positive integer. In Figs. 2 and 3 we show examples of ΩH with uniform and random distributions of ωs , s = 1, m, respectively. Random distributions of ωs , s = 1, m, are always subsets of underlying uniform distributions. We consider the case of square domain and square inclusions only for sake of simplicity. It is obvious that the proposed method is not limited to this case. In numerical tests we always choose min ds,t = 2H t ̸ =s
and
∫ x∈ωs y∈∂ Ω
|x − y| = H ,
i.e., minimal distance between neighboring subdomains ωs and ωt is not smaller than 2H and the minimal distance between ωs and ∂ Ω is not smaller than H, s, t = 1, m. The square meshes Ωh are designed by refinements of the meshes ΩH .
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
Fig. 3. An example of random distribution of ωs , s = 1, m, with d = 2H = Table 1 Accuracy of asymptotic expansions: d = 2H, H =
ε
10
−1
−2
10
10
1 , 128
h=
1 , 512
−3
1 , 16
h=H=
10−4
δu δ p(0)
6.3e−2 1.6e−1
8.0e−3 1.9e−2
8.0e−4 2.0e−3
8.3e−5 2.0e−4
δ u(1) δ p(1)
2.0e−3 4.3e−2
2.6e−4 5.4e−4
2.7e−6 5.6e−6
2.7e−8 5.6e−8
δ u(2) δ p(2)
6.9e−3 1.3e−4
9.0e−6 1.6e−5
9.2e−9 1.7e−8
9.3e−12 1.7e−11
δ u(3) δ p(3)
2.3e−3 4.1e−3
3.0e−7 5.3e−7
4.4e−11 7.5e−11
2.3e−13 2.9e−13
ε
10
−1
−2
10
10
1 , 256
m = 46.
m = 1024.
(0)
Table 2 Accuracy of asymptotic expansions: d = 2H, H =
1 , 32
9
h=
−3
1 , 512
m = 4096. 10−4
(0)
δu δ p(0)
7.3e−2 1.6e−1
9.6e−3 2.0e−2
9.0e−4 2.0e−3
9.8e−5 2.0e−4
δ u(1) δ p(1)
2.2e−2 4.4e−2
2.9e−4 5.4e−4
3.0e−6 5.5e−6
3.0e−8 5.5e−8
δ u(2) δ p(2)
7.0e−3 1.2e−2
9.1e−6 1.5e−5
9.4e−9 1.6e−8
9.2e−12 1.5e−11
δ u(3) δ p(3)
2.2e−3 3.7e−3
2.8e−7 4.6e−7
3.1e−11 5.0e−11
8.6e−14 4.0e−14
In numerical tests we compute the values
δu
(r)
u − u(r) A = , f −1 A
δp
(r)
p − p(r) B1 = . f −1
(36)
A
We solve the underlying systems (5), (22) with the relative accuracy 10−15 by the preconditioned Uzawa method, proposed and investigated in [4,5]. For solving the systems with the matrix Aε in (4) we use the PCG method with the algebraic multigrid preconditioner proposed in [7]. Concerning the results of iterative methods for system (5) we refer to [5]. 4.2. Convergence of asymptotic expansions In Tables 1–3 we show the accuracy results for uniform distributions of inclusions. Convergence condition (26) is satisfied, and numerical results clearly demonstrate good accuracy of the proposed expansion (19)–(20). We included
10
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 Table 3 Accuracy of asymptotic expansions: d = 2H, H =
ε
10
−1
10
−2
1 , 512
h=
1 , 512 −3
m = 16 384.
10
10−4
(0)
δu δ p(0)
9.3e−2 1.6e−1
8.0e−3 1.9e−2
1.0e−3 2.0e−3
1.0e−4 2.0e−4
δ u(1) δ p(1)
2.3e−2 4.0e−2
2.6e−4 5.4e−4
2.9e−6 5.0e−5
2.9e−8 5.0e−8
δ u(2) δ p(2)
6.6e−3 1.0e−2
9.5e−6 1.6-e5
8.3e−9 1.3e−8
8.4e−12 1.2e−11
δ u(3) δ p(3)
1.8e−3 2.6e−3
3.0e−7 5.3e−7
2.3e−11 3.4e−11
8.6e−14 5.4e−14
Fig. 4. Asymptotic distribution of ωs : ds = 14H, s = 1, m.
to the tables numerical results only for uniform distribution of ωs , s = 1, m, because the results for several underlying random distributions are absolutely the same. 4.3. Analysis of convergence estimates In this subsection, we investigate numerically accuracy of asymptotic expansions depending on the ratio of the distance between the neighboring inclusions ds,t to the size of ds . We are not able (too expensive) to compute the value of µm+1 in (25) but we are able to compute the minimal eigenvalue νs,min in (35) with given subdomains ω ˆ s , s = 1, m. In our case all the inclusions ωs have the same shape and size. To this end, we can compute the value of νs,min only for one particular inclusion ωs , 1 ⩽ s ⩽ m. We recall that due to (32) the estimate
µm+1 ⩾ min νs,min 1⩽s⩽m
holds. In numerical tests, we choose ω ˆ s as the square H-extension of ωs , i.e., ωˆ s are squares of the size dˆ s = In Fig. 4 we show a piece of ωs and ω ˆ s in the case k = 7. In Tables 4–6 we show accuracy results (36) depending on the values qε,h = ε/νmin . Recall that q=
k+2 ds , k
s = 1, m.
ε ε < qε,h = . µm+1 νmin
We can observe that accuracy deteriorates when qε,h becomes larger than one, and the expansions diverge for larger values of qε,h . 5. Conclusion and further remarks In this paper, we provided numerical analysis of the new asymptotic expansions for diffusion problems with highcontrast inclusions. These expansions were proposed and originally studied in [2]. The method and our analysis is not limited to the model problem described in Section 4.1. Numerical results in Section 4.2 are in very good agreement with theoretical estimates [2] briefly reviewed in Sections 2 and 3 of this paper. Also, numerical results in Tables 4–6 from Section 4.3 suggested that accuracy of asymptotic expansions deteriorates when distances between inclusions become smaller. In the future, we plan to extend our theoretical and numerical analysis to the 2D and 3D diffusion problems with complex shaped high-contrast inclusions, in particular, to problem with anisotropic inclusions.
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 Table 4 Accuracy of asymptotic expansions: ε = 10−1 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
126H
0.462
0.740
1.29
2.38
4.76
9.09
δ u(0) δ p(0)
8.7e−2 1.6e−1
5.6e−2 1.3e−1
4.0e−2 1.1e−1
3.3e−2 1.0e−1
3.8e−2 1.3e−1
4.0e−2 2.3e−1
δ p(1) δ u(1)
2.3e−2 4.0e−2
1.6e−2 3.0e−2
1.4e−2 2.9e−2
2.9e−2 6.9e−2
1.0e−1 3.3e−1
3.0e−1 1.35
δ u(2) δ p(2)
6.5e−3 1.0e−2
5.7e−3 9.4-e3
9.3e−3 1.7e−2
5.4e−2 1.2e−1
3.6e−1 1.15
2.06 8.99
δ u(3) δ p(3)
1.8e−3 2.6e−3
2.0e−3 3.2e−3
7.1e−3 1.4e−2
9.6e−2 2.3e−1
1.27 4.06
14.62 60.97
δ u(4) δ p(4)
5.0e−4 7.1e−4
7.6e−4 1.2e−3
6.6e−3 1.3e−2
1.5e−1 3.8e−1
3.96 12.95
107.25 424.47
qε,h =
ε νmin
Table 5 Accuracy of asymptotic expansions: ε = 10−2 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
126H
0.0462
00740
0.129
0.238
0.476
0.909
(0)
δu δ p(0)
1.0e−2 2.0e−2
7.0e−3 1.6e−2
5.0e−3 1.3e−2
4.7e−3 1.4e−2
9.6e−3 3.5e−2
2.0e−2 9.5e−2
δ p(1) δ u(1)
2.8e−4 4.9e−4
2.1e−4 3.8e−4
2.1e−4 4.2e−4
6.6e−4 1.5e−3
3.2e−3 1.1e−2
1.3e−2 6.2e−2
δ u(2) δ p(2)
8.1e−6 1.2e−5
7.4e−6 1.2e−5
1.5e−5 2.8e−5
1.1e−4 2.7e−4
1.1e−3 4.1e−3
9.7e−3 4.1e−2
δ u(3) δ p(3)
2.2e−7 3.3e−7
2.7e−7 4.3e−7
1.2e−6 2.5e−6
2.0e−5 5.0e−5
4.4e−4 1.4e−3
6.9e−3 2.8e−2
δ u(4) δ p(4)
6.2e−9 8.8e−9
1.0e−8 1.7e−8
1.9e−7 2.5e−7
3.5e−6 9.0e−6
1.3e−4 4.3e−4
5.2e−3 2.0e−2
qε,h =
ε
νmin
Table 6 Accuracy of asymptotic expansions: ε = 10−3 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
0.00462
0.00740
0.0129
0.0238
0.0476
0.0909
δ u(0) δ p(0)
1.0e−3 2.0e−3
7.0e−4 1.6e−3
5.1e−4 1.3e−3
5.3e−4 1.5e−3
1.3e−3 3.9e−3
3.2e−3 1.4e−2
δ p(1) δ u(1)
2.9e−6 5.0e−6
2.2e−6 3.9e−6
2.3e−6 4.5e−6
7.9e−6 1.9e−5
3.9e−5 1.4e−4
2.2e−4 9.6e−4
δ u(2) δ p(2)
8.3e−9 1.2e−8
7.7e−9 1.2e−8
1.6e−8 3.2e−8
1.3e−7 3.5e−7
1.3e−6 4.5e−6
1.5e−5 6.5e−5
δ u(3) δ p(3)
2.3e−11 3.2e−11
5.2e−11 7.3e−11
1.6e−10 3.0e−10
2.3e−9 5.8e−9
4.9e−8 1.7e−7
1.1e−6 4.5e−6
δ u(4) δ p(4)
2.5e−13 2.7e−13
3.1e−13 2.5e−13
1.2e−12 2.0e−12
4.1e−11 1.6e−10
1.8e−9 5.8e−9
8.8e−8 3.2e−7
qε,h =
ε νmin
11
126H
Acknowledgments This work was supported by the Russian Science Foundation through the grant 19-71-10094. References [1] V.M. Calo, Y. 465–494. [2] Y. Kuznetsov, [3] Y. Kuznetsov, [4] Y. Kuznetsov,
Efendiev, J. Galvis, Asymptotic expansions for high-contrast elliptic equations, Math. Model. Methods Appl. Sci. 24 (03) (2014) New homogenization method for diffusion equations, Russ. J. Numer. Anal. Math. Modell. 33 (2) (2018) 85–93. New iterative methods for singular perturbed positive definite matrices, Russ. J. Numer. Anal. Math. Modell. 15 (1) (2000) 65–71. Preconditioned iterative methods for algebraic saddle-point problems, J. Numer. Math. 17 (1) (2009) 67–75.
12
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
[5] Y. Gorb, V. Kramarenko, Y. Kuznetsov, Preconditioned iterative methods for diffusion problems with high-contrast inclusions, Numer. Linear Algebra Appl. 26 (4) (2019) e2243, http://dx.doi.org/10.1002/nla.2243. [6] O.B. Widlund, An Extension Theorem for Finite Element Spaces with Three Applications, in: Notes on Numerical Fluid Mechanics, vol. 16, Vieweg+Teubner Verlag, Wiesbaden, 1986, pp. 110–122. [7] Y. Kuznetsov, Multigrid domain decomposition method, Domain Decomposition Methods for PDEs, 1989, pp. 290–313.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
New asymptotic expansion for the diffusion problem with high-contrast inclusions: Numerical results ∗
V.K. Kramarenko a , , Yu. A. Kuznetsov b a b
Marchuk Institute of Numerical Mathematics, Russian Academy of Science, Gubkina 8, Moscow, Russian Federation University of Houston, Department of Mathematics, Houston, USA
article
info
Article history: Received 27 February 2019 Received in revised form 9 July 2019 Keywords: Diffusion equation High-contrast inclusions Finite elements Asymptotic expansions Saddle-point matrix
a b s t r a c t This paper is concerned with the new approach to the numerical solution of diffusion equations with high contrast coefficients. The classical P1 finite element system of equations is replaced by the equivalent system of the saddle point type. For the saddle point problem a new recently invented asymptotic expansion with respect to a small parameter is discussed and the estimates for the convergence factor are given. Numerical results completely confirm the theoretical accuracy and convergence results for the considered method. © 2019 Elsevier B.V. All rights reserved.
0. Introduction Elliptic problems with high contrast coefficients recently have been a topic of active research, see [1] and references therein. An approach based on the classical asymptotic expansion technique is among the most popular methods for solving this type of problems. A new approach to the asymptotic expansion for the P1 finite element approximations of the diffusion problem with high contrast inclusions was proposed and investigated in [2]. In what follows the main idea is to replace the classical P1 finite element problem by the equivalent problem of the saddle point type. This idea was invented originally for solving of algebraic systems with singularly perturbed symmetric positive definite matrices in [3]. In [4], this procedure was applied to solving P1 finite element discretizations of the diffusion equations with high contrast inclusions. An efficient preconditioned Lanczos method was proposed and investigated. Rigorous analysis with new convergence proof of the preconditioned iterative methods for underlying algebraic saddle point systems of equations was appeared in [5]. The current paper complements paper [2]. By virtue of the fact it contains the description of some practical aspects invented in [2] and numerical results which completely confirm the earlier theoretical results. The paper is organized as follows. In Section 1 we formulate the problem and discuss the algebraic aspects of the replacement of the original P1 finite element problem by the equivalent saddle type problem. In Section 2 we recall the new asymptotic expansion proposed in [2] and the easier convergence results. We discuss the estimates of the convergence factor for the new asymptotic expansion in Section 3. The model test problem and the numerical results are given in Section 4. For the sake of simplicity we consider the case of a square domain with multiple uniformly and randomly distributed square inclusions. The numerical results confirm the theoretical statements and the proposed estimates for the convergence factor. ∗ Corresponding author. E-mail addresses: [email protected] (V.K. Kramarenko), [email protected] (Yu. A. Kuznetsov). https://doi.org/10.1016/j.cam.2019.112363 0377-0427/© 2019 Elsevier B.V. All rights reserved.
2
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
1. Problem formulation We consider the diffusion problem
−∇ [k (x) ∇ u] = f , x ∈ Ω , u = 0, x ∈ ∂Ω, where Ω ∈ R2 is a polygonal domain with the boundary ∂ Ω and k (x) ⩾ 1 is a positive piecewise constant function. We introduce polygonal subdomains (inclusions) ωs of Ω , s = 1, m, where m is a positive integer, such that ωs ∩ωt = ∅ and ωs ∩ ∂ Ω = ∅, s, t = 1, m. For the sake of simplicity we assume that subdomains ωs are convex and
{ k (x) =
ks = 1 +
1,
1
εs
≡ const > 1 in ωs , s = 1, m,
otherwise,
and consider the weak formulation of : Find u ∈ H01 (Ω ) such that
∫ Ω
(k (x) ∇ u) · ∇v dx ≡
∫ ∇ u · ∇v dx + Ω
∫ ∫ m ∑ 1 ∇ u · ∇v dx = f v dx εs ωs Ω
(1)
s=1
∀u, v ∈ H01 (Ω ). Let Ωh be a triangular mesh in Ω conforming with the boundary ∂ωs , i.e. ∂ωs is the union of edges from Ωh , s = 1, m. We denote by Vh the classical P1-finite element subspace of H01 (Ω ) on Ωh . The P1-finite element method: Find uh ∈ Vh , such that
∫ Ω
(k (x) ∇ u) · ∇v =
∫
f v dx
Ω
∀v ∈ Vh
results in the algebraic system Aε u = f
(2)
with the N × N matrix Aε = Aε > 0 and the vector f ∈ R , where N is the number of interior mesh nodes in Ωh . With an appropriate ordering the matrix Aε can be presented as 2 × 2 block matrix T
[ Aε =
A11,ε A21
A12 A22
N
]
with n × n submatrix A11,ε = A11 + B1,ε , where B1,ε is the m × m block diagonal matrix with the diagonal blocks 1ε As ∈ Rns ×ns , s = 1, m. Here, ns is the number of ∑m nodes belonging to ωs , s = 1, m, n = s=1 ns , and matrices As are defined by the identities
(As us , v s ) =
∫ ωs
∇ uh · ∇vh dx ∀us , v s ∈ Rns ,
(3)
where uh , vh ∈ Vs,h and Vs,h is the restriction of Vh onto ωs , s = 1, m. In other words,
[ Aε =
]
A11 A21
[
A12 B + 1,ε A22 0
]
0 , 0
(4)
where
{ B1,ε = diag
1
ε1
A1 , . . . ,
1
εm
}
Am .
Using the procedure proposed in [3,4], we replace (2) by the system Au Bu
+ −
BT p Σε B1 p
= =
f, 0,
(5)
where
[
A11 A= A21
A12 A22
] (6)
is n × n matrix, B = [B1 0]
(7)
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
3
is n × N matrix, B1 = diag {A1 , . . . , Am } ,
Σε = diag {ε1 I1 , . . . , εm Im } , and p = Σε−1 0 u.
]
[
Here, Is are the identity ns × ns matrices, s = 1, m. It is obvious, that the matrix
(
A B
A=
)
BT
−Σε B1
is singular, dim(ker A) = m, ker BT = ker B1 , and any vector w ∈ ker B1 can be presented in the block form by
⎤ w1 ⎢ ⎥ w = ⎣ ... ⎦ , wm ⎡
where w s ∈ ker As are vectors with constant components, s = 1, m. The solution vector u ∈ Rn in (5) is unique, and the solution vector p ∈ Rn in (5) is unique up to an arbitrary additive vector w ∈ ker B1 . Henceforth, we shall always assume that the solution vector p in (5) is orthogonal to ker B1 , i.e. the uniqueness of p. Under the latter assumption p ⊥ ker B1 , system (5) is equivalent to the system
[ ] Lα
u ≡ p
(
BT −Σε B1 − Qα
A B
)[ ]
[ ]
u f , = p 0
(8)
where
{
}
Qα = diag α1 Q1 , . . . , αm Qm ,
αs = consts > 0, s = 1, m, Qs = Ms w s w Ts Ms ,
s = 1 , m,
and Ms are the underlying mass matrices defined by
(Ms us , v s ) =
∫ ωs
us,h vs,h dx
∀us,h , vs,h ∈ Vs,h ,
s = 1, m. We provide the following alternative definition of system (8). We replace the variational problem (1) by the equivalent variational saddle point problem: find u ∈ H01 (Ω ), ps ∈ H 1 (ωs ), s = 1, m, such that
∫ ∇ u · ∇v dx + Ω
∫ ωs
m ∫ ∑
ωs
s=1
∇ u · ∇ qs dx − εs
∫
∫ ωs
∇ ps · ∇v dx =
Ω
∇ ps · ∇ qs dx − αs
f v dx,
(9)
∫
∫ ωs
ps dx ·
ωs
qs dx = 0
(10)
∀v ∈ H01 (Ω ), qs ∈ H 1 (ωs ), s = 1, m, where αs are positive constants, s = 1, m.
Then system (8) results from the finite element discretization of (9)–(10): find uh ∈ Vh , ps,h ∈ Vs,h , s = 1, m, such that
∫ Ω
∫ ωs
∇ uh · ∇vh dx +
m ∫ ∑
ωs
s=1
∇ uh · ∇ qs,h dx − εs
∫
∫ ωs
∇ ps,h · ∇vh dx =
Ω
f vh dx,
∇ ps,h · ∇ qs,h dx − αs
∫
∫ ωs
ps,h dx ·
ωs
qs,h dx = 0
∀vh ∈ Vh , qs,h ∈ Vs,h , s = 1, m. 2. New asymptotic expansion 1 The matrix Lα in (8) is nonsingular for any εs ⩾ 0, and the entries of L− α are infinitely differentiable functions of εs ⩾ 0 1 and tend to the entries of L− as ε → 0, s = 1 , m. To this end, it is natural to use the classical asymptotic expansion for 0
4
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
the solution vectors u and p in (8) (the solution finite element functions uh and ps,h , s = 1, m, in (9)–(10)) with respect to parameters εs , s = 1, m. In this section, similar to [2], we consider only the case of ε1 = ε2 = · · · = εm = ε . First, we consider the system (8) with ε1 = ε2 = · · · = εm = 0, i.e. the system 0
Au 0 Bu
0
BT p 0 Qα p
+ −
f, 0,
= =
(11)
0
0
0
and estimate the errors u − u and p − p in A1 -norm and B1 -seminorm, respectively. We recall that Qα p = 0. Below we present some results and proofs from [2]. The error vectors satisfy the system
(
)(
BT
A B
−Σε B1
0
u−u 0 p−p
)
( =
0 Σε B1 p0
)
.
(12)
0
Eliminating the vector u − u in (12) we obtain the system
= −Σε B1 p0 ,
(S + Σε B1 ) p − p
) 0
(
(13)
where −1 S = BA−1 BT ≡ B1 S11 B1
and 1 S11 = A11 − A12 A− 22 A21 .
Consider the eigenvalue problem −1 B1 S11 B1 w = µB1 w
(14)
with the eigenvalues 0 = µ1 = µ2 = · · · = µm < µm+1 ⩽ 1 where µm+1 is the minimal positive eigenvalue. The inequality µm+1 ⩽ 1 was proved in [3]. Then from (13) we derive the estimate
p − p0
B1
⩽
ε
· p0 B ,
µm+1
(15)
1
where ∥·∥B1 denotes 0 the semi-norm generated by the matrix B1 . 0 To estimate p B in (15), we return to system (11). Eliminating the vector u we have the equation 1
0
−1
Sp = −B1 A
f.
Then the estimate
0 p
B1
=
1
µm+1
f
(16)
A− 1
follows from the fact that
1 12 − 1 B A 2 = A− 21 B 2 ⩽ 1 1 1 2
(17)
2
where ∥·∥A−1 is the norm generated by A−1 . Thus, we obtain the estimate
p − p0
B1
⩽
ε · f A−1 . µ2m+1 0
Multiplying the first equation in (12) by the vector u − u , and applying the Schwarz inequality and (17), we get the estimates
u − u0 ⩽ p − p 0 ⩽ A B 1
ε µ
2 m+1
· f A− 1 .
(18)
We now consider the expansions u=
∞ ∑
ε i u(i) ,
(19)
ε i p(i) ,
(20)
i=0
p=
∞ ∑ i=0
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 (i)
where the vectors u , p (i)
Au (i) Bu
+
T
B p
(i)
(i)
5
satisfy Eq. (11) for i = 0 and the equations
= =
0, (i−1) B1 p ,
i = 1, 2, . . .. Similar to above we readily derive the estimates
(i) p
B1
⩽
1
µm+1
· p(i−1) B ⩽ 1
1
(µm+1 )i
· p(0) B , 1
(21)
i = 1, 2, . . .. Define the error vectors (r)
y
=u−
r ∑
ε l u(l) ,
l=0
z
(r)
=p−
r ∑
ε l p(l)
l=0
for some integer r ⩾ 1. (r) (r) It is obvious that the vectors y and z satisfy the system (r)
Ay (r) By
+ −
(r)
BT z ε B1 z (r)
= =
0, −εr +1 B1 p(r) ,
(22)
r = 1, 2, . . .. (r) Again eliminating the vector y in (22) we obtain the system
(S + ε B1 ) z (r) = −ε r +1 B1 p(r) . Now, using the estimate
(r) z
B1
⩽
ε r +1 · p(r) B , 1 µm+1
and (16) and (21), and the estimate
(r) y ⩽ z (r) , A B
(23)
1
we obtain the final estimates r +1 r y ⩽ z (r) ⩽ ε · f A−1 , r +2 A B1 µm+1
r = 0, 1, . . . .
(24)
On the bases of estimates (18) and (24) we reach the following conclusion: Statement 1. The estimate r +1 (r) y ⩽ ε · u∗ A r +2 A µm+1
(25)
∗
holds for any r ⩾ 0 where u = A−1 f . Here, u∗h is the P1-finite element solution of the Poisson problem
− △ u∗ = f in Ω , u∗ = 0 on ∂ Ω on the mesh Ωh . Statement 2. The condition
ε
µm+1
<1
(26) (r)
is sufficient for convergence of uh to uh as r −→ +∞. 3. Estimation from below for µm+1 We replace (14) by an equivalent eigenvalue problem: Find µ ∈ R, v ∈ RN , w ∈ Rn , such that Av Bv
+
BT w
= =
0, −µB1 w,
(27)
6
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
where
v = −A−1 BT w ≡
[ ] v1 v2
−1 and v 1 = −S11 B1 w . It follows that
A21 v 1 + A22 v 2 = 0, i.e., v2,h ∈ VΩ \ω,h is the h-harmonic extension of v1,h ∈ Vω,h . Here Vω,h and VΩ \ω,h are the restrictions of Vh onto ω and Ω \ ω, respectively. The finite element formulation of (27) is as follows: Find µ ∈ R, vh ∈ Vh and wh ∈ Vω,h , such that
∫
∫
∫Ω ω
∇vh · ∇ϕh dx
+ ω
∇wh · ∇ϕh dx
0,
=
(28)
∫ ∇vh · ∇ψh dx
=
−µ ω
∇wh · ∇ψh dx
for any ϕh ∈ Vh , ψh ∈ Vω,h . To consider only nonzero eigenvalues, we impose additional conditions
∫
wh dx = 0,
ωs
s = 1, m.
(29)
It is clear that µm+1 in (14) is the minimal eigenvalue in (28)–(29). It was established in [4,5] with the reference to the finite element norm preserving extension result in [6] that in the case of regular shaped quasiuniform meshes Ωh the value of µm+1 is bounded form below by a positive constant independent of the mesh. Moreover, it was further established that under the additional assumptions that ωs are convex and regular shaped, and ds,t = inf |x − y| ⩾ αs ds ,
(30)
x∈ωs , y∈ωt , t ̸ =s
where αs are the positive constants and ds = |ωs |1/2 , s, t = 1, m, the value of µm+1 is bounded from below by a positive constant also independent of ds , s = 1, m. An alternative proof of the above results was recently proposed in [2]. In this section we follow the latest proof. Let µ be an eigenvalue and vh , wh corresponding eigenfunctions in (28)–(29). If we choose ϕh = vh and ψh = wh , we obtain the equality
∫
|∇vh |2 dx = µ
Ω
∫ ω
|∇wh |2 dx.
Taking ϕh = vh and using the inequality
⏐ [∫ ⏐∫ ]1/2 [∫ ]1/2 ⏐ ⏐ 2 2 ⏐ ∇wh · ∇vh dx⏐ ⩽ |∇w | |∇v | dx dx h h ⏐ ⏐ ω
ω
Ω
we obtain
∫ Ω
|∇vh |2 dx ⩽
∫ ω
|∇wh |2 dx
and conclude that the eigenvalues µ in (28)–(29) belong to the segment (0; 1]. Now we choose ϕh ∈ Vh such that ϕh = wh in ω and ϕh = w2,h in Ω \ ω, where w2,h is a finite element extension of wh from ω into Ω \ ω, i.e. w2,h = wh on ∂ω and w2,h ∈ VΩ \ω,h . From the equality
∫
∫ Ω
∇vh · ∇ϕh dx = −
ω
|∇wh |2 dx,
we derive the estimate
∫ ω
[ ] |∇wh |2 dx ⩽ 1 + C 2
∫ Ω
|∇vh |2 dx
with
∫ C2 =
Ω \ω
∫ ω
⏐ ⏐ ⏐∇w2,h ⏐2 dx . |∇wh |2 dx
(31)
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
7
Fig. 1. An example of ωs and ω ˆ s.
The value of C 2 is minimal when w2,h is the h-harmonic finite element extension of wh from ω to Ω \ ω. Choosing φh = vh and ψh = wh in (28) and using (31) we obtain the estimate 1
µm+1
⩽ 1 + C 2.
To estimate C 2 from above, we construct an appropriate function w2,h . We define in Ω a set of nonoverlapping open polygons ω ˆ s , s = 1, m, with the boundaries ∂ ωˆ s , such that ωs ∈ ωˆ s , i.e., ∂ωs ∩ ∂ ωˆ s = 0 and ∂ ωˆ s ∩ ∂ Ω = 0, s = 1, m. We assume that ω ˆ s are conforming with the mesh Ωh , i.e., ∂ ωˆ s are the unions of sides of triangular mesh cells in Ωh , s = 1, m. An example of hexagonal ωs and ω ˆ s , 1 ⩽ s ⩽ m, is given in Fig. 1, where the boundary ∂ ωˆ s is denoted by dash lines. We subject w2,h to the conditions
w2,h = 0 in Ω \ ωˆ ⋃ where ω ˆ = m ˆ s, s=1 ω w2,h = wh on ∂ωs , s = 1, m (s)
and w2,h ≡ w2,h in ω ˆ s , where
∫ ωˆ s
∇w2(s),h · ∇vh dx = 0 ∀vh ∈ Vωˆ s \ωs ,
(s) 2,h
ˆ s , s = 1, m. i.e., w are the h-harmonic extensions of wh from ωs into ω We additionally assume that ∂ ω ˆ s satisfy the conditions β · ds ⩽ inf |x − y| x∈∂ωs y∈∂ ω ˆs
with positive constants β independent of Ωh , and the meshes ω ˆ s \ ωs are regular shaped and quasiuniform, s = 1, m. Under these assumptions, it was shown in [2,4,5] that the values of 2
∫ cs2
=
ωˆ s \ωs
∫
ωs
|∇w2(s),h | dx |∇wh |2 dx
,
s = 1, m, are bounded from above by positive constants independent of the mesh Ωh and values of ds , s = 1, m. It follows immediately that 1
µm+1
⩽ 1 + max cs2 .
(32)
1⩽s⩽m
The values of 1/(1 + cs2 ) are the minimal eigenvalues of the eigenvalue problems: Find νs ∈ R, vs,h ∈ Vωˆ s ,h , ws,h ∈ Vωs ,h , such that
∫ ∫ ωˆ s ∇vs,h · ∇φh dx ∇vs,h · ∇ψh dx ωˆ s
+
∫
ωs
∇ws,h · ∇φh dx
= =
−νs
0, ∇w s,h · ∇ψs,h dx ω
∫
(33)
8
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
Fig. 2. An example of uniform distribution of ωs , s = 1, m, with d = 2H =
1 , 8
h=H=
1 , 16
m = 16.
∀φh ∈ Vωˆ s ,h , ψh ∈ Vωs ,h under the condition ∫ ws,h dx = 0,
(34)
ωs
s = 1, m. The eigenvalue problem (33), (34) is equivalent to the algebraic eigenvalue problem A(s) us B(s) v s
+
[
B(s)
]T
ws
= =
0, −νs As v s ,
(35)
under the condition (34), where As are defined in (3) and A(s) and B(s) are the corresponding blocks of the matrices A and B in (6) and (7), respectively, s = 1, m. For instance, B(s) = As
[
0 ,
]
s = 1, m. The minimal eigenvalues in (33)–(35) can be computed numerically, or estimated, s = 1, m. 4. Numerical results 4.1. Model test problem Let Ω be the unit square and ΩH be the square mesh with the mesh size H =
1 NH
, where NH is a positive integer. We
define the set of ωs , s = 1, m, assuming that each ωs is a square and a union of mesh cells in ΩH , 1 ⩽ s ⩽ m. In numerical tests we always assume that ωs have the same size d ≡ ds = kH, s = 1, m, where k ⩾ 2 is a positive integer. In Figs. 2 and 3 we show examples of ΩH with uniform and random distributions of ωs , s = 1, m, respectively. Random distributions of ωs , s = 1, m, are always subsets of underlying uniform distributions. We consider the case of square domain and square inclusions only for sake of simplicity. It is obvious that the proposed method is not limited to this case. In numerical tests we always choose min ds,t = 2H t ̸ =s
and
∫ x∈ωs y∈∂ Ω
|x − y| = H ,
i.e., minimal distance between neighboring subdomains ωs and ωt is not smaller than 2H and the minimal distance between ωs and ∂ Ω is not smaller than H, s, t = 1, m. The square meshes Ωh are designed by refinements of the meshes ΩH .
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
Fig. 3. An example of random distribution of ωs , s = 1, m, with d = 2H = Table 1 Accuracy of asymptotic expansions: d = 2H, H =
ε
10
−1
−2
10
10
1 , 128
h=
1 , 512
−3
1 , 16
h=H=
10−4
δu δ p(0)
6.3e−2 1.6e−1
8.0e−3 1.9e−2
8.0e−4 2.0e−3
8.3e−5 2.0e−4
δ u(1) δ p(1)
2.0e−3 4.3e−2
2.6e−4 5.4e−4
2.7e−6 5.6e−6
2.7e−8 5.6e−8
δ u(2) δ p(2)
6.9e−3 1.3e−4
9.0e−6 1.6e−5
9.2e−9 1.7e−8
9.3e−12 1.7e−11
δ u(3) δ p(3)
2.3e−3 4.1e−3
3.0e−7 5.3e−7
4.4e−11 7.5e−11
2.3e−13 2.9e−13
ε
10
−1
−2
10
10
1 , 256
m = 46.
m = 1024.
(0)
Table 2 Accuracy of asymptotic expansions: d = 2H, H =
1 , 32
9
h=
−3
1 , 512
m = 4096. 10−4
(0)
δu δ p(0)
7.3e−2 1.6e−1
9.6e−3 2.0e−2
9.0e−4 2.0e−3
9.8e−5 2.0e−4
δ u(1) δ p(1)
2.2e−2 4.4e−2
2.9e−4 5.4e−4
3.0e−6 5.5e−6
3.0e−8 5.5e−8
δ u(2) δ p(2)
7.0e−3 1.2e−2
9.1e−6 1.5e−5
9.4e−9 1.6e−8
9.2e−12 1.5e−11
δ u(3) δ p(3)
2.2e−3 3.7e−3
2.8e−7 4.6e−7
3.1e−11 5.0e−11
8.6e−14 4.0e−14
In numerical tests we compute the values
δu
(r)
u − u(r) A = , f −1 A
δp
(r)
p − p(r) B1 = . f −1
(36)
A
We solve the underlying systems (5), (22) with the relative accuracy 10−15 by the preconditioned Uzawa method, proposed and investigated in [4,5]. For solving the systems with the matrix Aε in (4) we use the PCG method with the algebraic multigrid preconditioner proposed in [7]. Concerning the results of iterative methods for system (5) we refer to [5]. 4.2. Convergence of asymptotic expansions In Tables 1–3 we show the accuracy results for uniform distributions of inclusions. Convergence condition (26) is satisfied, and numerical results clearly demonstrate good accuracy of the proposed expansion (19)–(20). We included
10
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 Table 3 Accuracy of asymptotic expansions: d = 2H, H =
ε
10
−1
10
−2
1 , 512
h=
1 , 512 −3
m = 16 384.
10
10−4
(0)
δu δ p(0)
9.3e−2 1.6e−1
8.0e−3 1.9e−2
1.0e−3 2.0e−3
1.0e−4 2.0e−4
δ u(1) δ p(1)
2.3e−2 4.0e−2
2.6e−4 5.4e−4
2.9e−6 5.0e−5
2.9e−8 5.0e−8
δ u(2) δ p(2)
6.6e−3 1.0e−2
9.5e−6 1.6-e5
8.3e−9 1.3e−8
8.4e−12 1.2e−11
δ u(3) δ p(3)
1.8e−3 2.6e−3
3.0e−7 5.3e−7
2.3e−11 3.4e−11
8.6e−14 5.4e−14
Fig. 4. Asymptotic distribution of ωs : ds = 14H, s = 1, m.
to the tables numerical results only for uniform distribution of ωs , s = 1, m, because the results for several underlying random distributions are absolutely the same. 4.3. Analysis of convergence estimates In this subsection, we investigate numerically accuracy of asymptotic expansions depending on the ratio of the distance between the neighboring inclusions ds,t to the size of ds . We are not able (too expensive) to compute the value of µm+1 in (25) but we are able to compute the minimal eigenvalue νs,min in (35) with given subdomains ω ˆ s , s = 1, m. In our case all the inclusions ωs have the same shape and size. To this end, we can compute the value of νs,min only for one particular inclusion ωs , 1 ⩽ s ⩽ m. We recall that due to (32) the estimate
µm+1 ⩾ min νs,min 1⩽s⩽m
holds. In numerical tests, we choose ω ˆ s as the square H-extension of ωs , i.e., ωˆ s are squares of the size dˆ s = In Fig. 4 we show a piece of ωs and ω ˆ s in the case k = 7. In Tables 4–6 we show accuracy results (36) depending on the values qε,h = ε/νmin . Recall that q=
k+2 ds , k
s = 1, m.
ε ε < qε,h = . µm+1 νmin
We can observe that accuracy deteriorates when qε,h becomes larger than one, and the expansions diverge for larger values of qε,h . 5. Conclusion and further remarks In this paper, we provided numerical analysis of the new asymptotic expansions for diffusion problems with highcontrast inclusions. These expansions were proposed and originally studied in [2]. The method and our analysis is not limited to the model problem described in Section 4.1. Numerical results in Section 4.2 are in very good agreement with theoretical estimates [2] briefly reviewed in Sections 2 and 3 of this paper. Also, numerical results in Tables 4–6 from Section 4.3 suggested that accuracy of asymptotic expansions deteriorates when distances between inclusions become smaller. In the future, we plan to extend our theoretical and numerical analysis to the 2D and 3D diffusion problems with complex shaped high-contrast inclusions, in particular, to problem with anisotropic inclusions.
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363 Table 4 Accuracy of asymptotic expansions: ε = 10−1 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
126H
0.462
0.740
1.29
2.38
4.76
9.09
δ u(0) δ p(0)
8.7e−2 1.6e−1
5.6e−2 1.3e−1
4.0e−2 1.1e−1
3.3e−2 1.0e−1
3.8e−2 1.3e−1
4.0e−2 2.3e−1
δ p(1) δ u(1)
2.3e−2 4.0e−2
1.6e−2 3.0e−2
1.4e−2 2.9e−2
2.9e−2 6.9e−2
1.0e−1 3.3e−1
3.0e−1 1.35
δ u(2) δ p(2)
6.5e−3 1.0e−2
5.7e−3 9.4-e3
9.3e−3 1.7e−2
5.4e−2 1.2e−1
3.6e−1 1.15
2.06 8.99
δ u(3) δ p(3)
1.8e−3 2.6e−3
2.0e−3 3.2e−3
7.1e−3 1.4e−2
9.6e−2 2.3e−1
1.27 4.06
14.62 60.97
δ u(4) δ p(4)
5.0e−4 7.1e−4
7.6e−4 1.2e−3
6.6e−3 1.3e−2
1.5e−1 3.8e−1
3.96 12.95
107.25 424.47
qε,h =
ε νmin
Table 5 Accuracy of asymptotic expansions: ε = 10−2 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
126H
0.0462
00740
0.129
0.238
0.476
0.909
(0)
δu δ p(0)
1.0e−2 2.0e−2
7.0e−3 1.6e−2
5.0e−3 1.3e−2
4.7e−3 1.4e−2
9.6e−3 3.5e−2
2.0e−2 9.5e−2
δ p(1) δ u(1)
2.8e−4 4.9e−4
2.1e−4 3.8e−4
2.1e−4 4.2e−4
6.6e−4 1.5e−3
3.2e−3 1.1e−2
1.3e−2 6.2e−2
δ u(2) δ p(2)
8.1e−6 1.2e−5
7.4e−6 1.2e−5
1.5e−5 2.8e−5
1.1e−4 2.7e−4
1.1e−3 4.1e−3
9.7e−3 4.1e−2
δ u(3) δ p(3)
2.2e−7 3.3e−7
2.7e−7 4.3e−7
1.2e−6 2.5e−6
2.0e−5 5.0e−5
4.4e−4 1.4e−3
6.9e−3 2.8e−2
δ u(4) δ p(4)
6.2e−9 8.8e−9
1.0e−8 1.7e−8
1.9e−7 2.5e−7
3.5e−6 9.0e−6
1.3e−4 4.3e−4
5.2e−3 2.0e−2
qε,h =
ε
νmin
Table 6 Accuracy of asymptotic expansions: ε = 10−3 , H = ds
1 , 512
h=
1 . 512
2H
6H
14H
30H
62H
0.00462
0.00740
0.0129
0.0238
0.0476
0.0909
δ u(0) δ p(0)
1.0e−3 2.0e−3
7.0e−4 1.6e−3
5.1e−4 1.3e−3
5.3e−4 1.5e−3
1.3e−3 3.9e−3
3.2e−3 1.4e−2
δ p(1) δ u(1)
2.9e−6 5.0e−6
2.2e−6 3.9e−6
2.3e−6 4.5e−6
7.9e−6 1.9e−5
3.9e−5 1.4e−4
2.2e−4 9.6e−4
δ u(2) δ p(2)
8.3e−9 1.2e−8
7.7e−9 1.2e−8
1.6e−8 3.2e−8
1.3e−7 3.5e−7
1.3e−6 4.5e−6
1.5e−5 6.5e−5
δ u(3) δ p(3)
2.3e−11 3.2e−11
5.2e−11 7.3e−11
1.6e−10 3.0e−10
2.3e−9 5.8e−9
4.9e−8 1.7e−7
1.1e−6 4.5e−6
δ u(4) δ p(4)
2.5e−13 2.7e−13
3.1e−13 2.5e−13
1.2e−12 2.0e−12
4.1e−11 1.6e−10
1.8e−9 5.8e−9
8.8e−8 3.2e−7
qε,h =
ε νmin
11
126H
Acknowledgments This work was supported by the Russian Science Foundation through the grant 19-71-10094. References [1] V.M. Calo, Y. 465–494. [2] Y. Kuznetsov, [3] Y. Kuznetsov, [4] Y. Kuznetsov,
Efendiev, J. Galvis, Asymptotic expansions for high-contrast elliptic equations, Math. Model. Methods Appl. Sci. 24 (03) (2014) New homogenization method for diffusion equations, Russ. J. Numer. Anal. Math. Modell. 33 (2) (2018) 85–93. New iterative methods for singular perturbed positive definite matrices, Russ. J. Numer. Anal. Math. Modell. 15 (1) (2000) 65–71. Preconditioned iterative methods for algebraic saddle-point problems, J. Numer. Math. 17 (1) (2009) 67–75.
12
V.K. Kramarenko and Yu. A. Kuznetsov / Journal of Computational and Applied Mathematics 365 (2020) 112363
[5] Y. Gorb, V. Kramarenko, Y. Kuznetsov, Preconditioned iterative methods for diffusion problems with high-contrast inclusions, Numer. Linear Algebra Appl. 26 (4) (2019) e2243, http://dx.doi.org/10.1002/nla.2243. [6] O.B. Widlund, An Extension Theorem for Finite Element Spaces with Three Applications, in: Notes on Numerical Fluid Mechanics, vol. 16, Vieweg+Teubner Verlag, Wiesbaden, 1986, pp. 110–122. [7] Y. Kuznetsov, Multigrid domain decomposition method, Domain Decomposition Methods for PDEs, 1989, pp. 290–313.