Local and Parallel Finite Element Algorithm Based On Multilevel Discretization for Eigenvalue Problem

0.26 MB·

by  Yu Li

Local and Parallel Finite Element Algorithm Based
On Multilevel Discretization for Eigenvalue
Problem∗
4
1 Yu Li†, Xiaole Han‡, Hehu Xie§ and Chunguang You¶
0
2

n
a
J
0 Abstract
2
A local and parallel algorithm based on the multilevel discretization is

]
proposed in this paper to solve the eigenvalue problem by the finite element
A
method. With this new scheme, solving the eigenvalue problem in the finest
N
gridistransferredtosolutionsoftheeigenvalueproblemsonthecoarsestmesh
.
h
and a series of solutions of boundary value problems by using the local and
t
a parallel algorithm. The computational work in each processor can reach the
m
optimalorder. Therefore,thistypeofmultilevel localandparallelmethodim-
[
provestheoverallefficiencyofsolvingtheeigenvalueproblem. Somenumerical

1 experiments are presented to validate the efficiency of the new method.
v
9 Keywords. eigenvalue problem, multigrid, multilevel correction, local
6
and parallel method, finite element method.
9
4
AMS subject classifications. 65N30, 65N25, 65L15, 65B99.
.
1
0
4
1 1 Introduction
:
v
i
X Solving large scale eigenvalue problems becomes a fundamental problem in modern
r science and engineering society. However, it is always a very difficult task to solve
a
∗This work was supported in part by National Science Foundations of China (NSFC 11001259,
11371026, 11201501, 11031006, 2011CB309703), the National Center for Mathematics and Inter-
disciplinary Science, CAS and the President Foundation of AMSS-CAS.
†LSEC,ICMSEC,AcademyofMathematicsandSystemsScience,ChineseAcademyofSciences,
Beijing 100190,China ([email protected])
‡LSEC,ICMSEC,AcademyofMathematicsandSystemsScience,ChineseAcademyofSciences,
Beijing 100190,China ([email protected])
§LSEC,ICMSEC,AcademyofMathematicsandSystemsScience,ChineseAcademyofSciences,
Beijing 100190,China ([email protected])
¶LSEC,ICMSEC,AcademyofMathematicsandSystemsScience,ChineseAcademyofSciences,
Beijing 100190,China ([email protected])
1
high-dimensional eigenvalue problems which come from physical and chemistry sci-
ences. Xu and Zhou [23] give a type of two-grid discretization method to improve
the efficiency of the solution of eigenvalue problems. By the two-grid method, the
solution of eigenvalue problem on a fine mesh is reduced to a solution of eigenvalue
problem on a coarse mesh (depends on the fine mesh) and a solution of the corre-
sponding boundary value problem on the fine mesh [23]. For more details, please
read [20, 21]. Combing the two-grid idea and the local and parallel finite element
technique [22], a type of local and parallel finite element technique to solve the
eigenvalue problems is given in [24] (also see [9]). Recently, a new type of multilevel
correction method for solving eigenvalue problems which can be implemented on
multilevel grids is proposed in [12]. In the multilevel correction scheme, the solution
of eigenvalue problem on a finest mesh can be reduced to a series of solutions of
the eigenvalue problem on a very coarse mesh (independent of finest mesh) and a
series of solutions of the boundary value problems on the multilevel meshes. The
multilevel correction method gives a way to construct a type of multigrid scheme
for the eigenvalue problem [13].
Inthispaper, weproposeatypeofmultilevel localandparallel scheme tosolve the
eigenvalue problem based on the combination of the multilevel correction method
and the local and parallel technique. An special property of this scheme is that we
can do the local and parallel computing for any times and then the mesh size of
original coarse triangulation is independent of the finest triangulation. With this
new method, the solution of the eigenvalue problem will not be more difficult than
the solution of the boundary value problems by the local and parallel algorithm
since the main part of the computation in the multilevel local and parallel method
is solving the boundary value problems.
The standard Galerkin finite element method for eigenvalue problem has been
extensively investigated, e.g. Babuˇska and Osborn [2, 3], Chatelin [7] and references
cited therein. There also exists analysis for the local and parallel finite element
method for the boundary value problems and eigenvalue problems [9, 16, 17, 22,
23, 24]. Here we adopt some basic results in these papers for our analysis. The
corresponding error and computational work estimates of the proposed multilevel
local and parallel scheme for the eigenvalue problem will be analyzed. Based on the
analysis, the new method can obtain optimal errors with an optimal computational
work in each processor.
An outline of this paper goes as follows. In the next section, a basic theory about
the local error estimate of the finite element method is introduced. In Section 3, we
introduce the finite element method for the eigenvalue problem and the correspond-
ing error estimates. A local and parallel type of one correction step and multilevel
correction algorithm will be given in Section 4. The estimate of the computational
work for the multilevel local and parallel algorithm is presented in section 5. In
Section 6, two numerical examples are presented to validate our theoretical analysis
and some concluding remarks are given in the last section.
2
2 Discretization by finite element method
In this section, we introduce some notation and error estimates of the finite element
approximation for linear elliptic problem. The letter C (with or without subscripts)
denotesagenericpositiveconstantwhichmaybedifferentatitsdifferentoccurrences
through the paper. For convenience, the symbols ., & and ≈ will be used in this
paper. That x . y ,x & y and x ≈ y , mean that x ≤ C y , x ≥ c y and
1 1 2 2 3 3 1 1 1 2 2 2
c x ≤ y ≤ C x for some constants C ,c ,c and C that are independent of mesh
3 3 3 3 3 1 2 3 3
sizes (see, e.g., [19]). We shall use the standard notation for Sobolev spaces Ws,p(Ω)
and their associated norms and seminorms (see, e.g., [1]). For p = 2, we denote
Hs(Ω) = Ws,2(Ω) and H1(Ω) = {v ∈ H1(Ω) : v| = 0}, where v| = 0 is in the
0 ∂Ω ∂Ω
sense of trace, k·k = k·k .
s,Ω s,2,Ω
For G ⊂ D ⊂ Ω, the notation G ⊂⊂ D means that dist(∂D \∂Ω,∂G \∂Ω) > 0
(see Figure 1). It is well known that any w ∈ H1(Ω ) can be naturally extended to
0 0
be a function in H1(Ω) with zero outside of Ω , where Ω ⊂ Ω. Thus we will show
0 0 0
this fact by the abused notation H1(Ω ) ⊂ H1(Ω).
0 0 0
Figure 1: G ⊂⊂ D ⊂⊂ Ω
2.1 Finite element space
Now, let us define the finite element space. First we generate a shape-regular de-
composition T (Ω) of the computing domain Ω ⊂ Rd (d = 2,3) into triangles or
h
rectangles for d = 2 (tetrahedrons or hexahedrons for d = 3). The diameter of a
cell K ∈ T (Ω) is denoted by h . The mesh size function is denoted by h(x) whose
h K
value is the diameter h of the element K including x.
K
For generality, following [22, 24], we shall consider a class of finite element spaces
that satisfy certain assumptions. Now we describe such assumptions.
A.0. There exists γ > 1 such that
hγ . h(x), ∀x ∈ Ω,
Ω
where h = max h(x) is the largest mesh size of T (Ω).
Ω x∈Ω h
3
Based on the triangulation T (Ω), we define the finite element space V (Ω) as
h h
follows
V (Ω) = v ∈ C(Ω¯) : v| ∈ P , ∀K ∈ T (Ω) ,
h K k h
(cid:8) (cid:9)
where P denotes the space of polynomials of degree not greater than a positive
k
integer k. Then we know V (Ω) ⊂ H1(Ω) and define V (Ω) = V (Ω) ∩ H1(Ω).
h 0h h 0
Given G ⊂ Ω, we define V (G) and T (G) to be the restriction of V (Ω) and T (Ω)
h h h h
to G, respectively, and
V (G) = v ∈ V (Ω) : suppv ⊂⊂ G .
0h h
(cid:8) (cid:9)
For any G ⊂ Ω mentioned in this paper, we assume that it aligns with the partition
T (Ω).
h
As we know, the finite element space V satisfy the following proposition (see,
h
e.g., [6, 8, 22, 24]).
Proposition 2.1. (Fractional Norm) For any G ⊂ Ω, we have
inf kw−vk . kwk , ∀w ∈ V (Ω). (2.1)
1,G 1/2,∂G h
v∈V0h(G)
2.2 A linear elliptic problem
In this subsection, we repeat some basic properties of a second order elliptic bound-
ary value problem and its finite element discretization, which will be used in this
paper. The following results is presented in [16, 17, 22, 24].
We consider the homogeneous boundary value problem
Lu = f, in Ω,
(2.2)
(cid:26) u = 0, on ∂Ω.
Here the linear second order elliptic operator L : H1(Ω) → H−1(Ω) is define as
0
Lu = −div(A∇u),
where A = (a ) ∈ Rd×d is uniformly positive definite symmetric on Ω with
ij 1≤i,j≤d
a ∈ W1,∞(Ω). The weak form for (2.2) is as follows:
ij
Find u ≡ L−1f ∈ H1(Ω) such that
0
a(u,v) = (f,v), ∀v ∈ H1(Ω), (2.3)
0
where (·,·) is the standard inner-product of L2(Ω) and
a(u,v) = A∇u,∇v .
(cid:0) (cid:1)
4
As we know
kwk2 . a(w,w), ∀w ∈ H1(Ω).
1,Ω 0
We assume (c.f. [11]) that the following regularity estimate holds for the solution of
(2.2) or (2.3)
kuk . kfk
1+α,Ω −1+α,Ω
for some α ∈ (0,1] depending on Ω and the coefficient of L.
For some G ⊂ Ω, we need the following regularity assumption
R(G). For any f ∈ L2(G), there exists a w ∈ H1(G) satisfying
0
a(v,w) = (f,v), ∀v ∈ H1(G)
0
and
kuk . kfk .
1+α,G −1+α,G
Fortheanalysis, wedefinetheGalerkin-ProjectionoperatorP : H1(Ω) → V (Ω)
h 0 0h
by
a(u−P u,v) = 0, ∀v ∈ V (Ω) (2.4)
h 0h
and apparently
kP uk . kuk , ∀u ∈ H1(Ω). (2.5)
h 1,Ω 1,Ω 0
Basedon(2.5),theglobalpriorierrorestimatecanbeobtainedfromtheapproximate
properties of the finite dimensional subspace V (Ω) (cf. [6, 8]). For the following
0h
analysis, we introduce the following quantity:
ρ (h) = sup inf kL−1f −vk . (2.6)
Ω 1,Ω
f∈L2(Ω),kfk0,Ω=1v∈V0h(Ω)
Similarly, we can also define ρ (h) if Assumption R(G) holds.
G
The following results can be found in [3, 6, 8, 23, 24].
Proposition 2.2.
k(I −P )L−1fk . ρ (h)kfk , ∀f ∈ L2(Ω),
h 1,Ω Ω 0,Ω
ku−P uk . ρ (h)ku−P uk , ∀u ∈ H1(Ω).
h 0,Ω Ω h 1,Ω 0
Now, we state an important and useful result about the local error estimates
[16, 17, 24] which will be used in the following.
Proposition 2.3. Suppose that f ∈ H−1(Ω) and G ⊂⊂ Ω ⊂ Ω. If Assumptions
0
A.0 holds and w ∈ V (Ω ) satisfies
h 0
a(w,v) = (f,v), ∀v ∈ V (Ω ).
0h 0
Then we have the following estimate
kwk . kwk +kfk .
1,G 0,Ω0 −1,Ω0
5
3 Error estimates for eigenvalue problems
In this section, we introduce the concerned eigenvalue problem and the correspond-
ing finite element discretization.
In this paper, we consider the following eigenvalue problem:
Find (λ,u) ∈ R×H1(Ω) such that b(u,u) = 1 and
0
a(u,v) = λb(u,v), ∀v ∈ H1(Ω), (3.1)
0
where
b(u,u) = (u,u).
For the eigenvalue λ, there exists the following Rayleigh quotient expression (see,
e.g., [2, 3, 23])
a(u,u)
λ = .
b(u,u)
From [3, 7], we know the eigenvalue problem (3.1) has an eigenvalue sequence {λ } :
j
0 < λ ≤ λ ≤ ··· ≤ λ ≤ ··· , lim λ = ∞,
1 2 k k
k→∞
and the associated eigenfunctions
u ,u ,··· ,u ,··· ,
1 2 k
where b(u ,u ) = δ . In the sequence {λ }, the λ are repeated according to their
i j ij j j
multiplicity.
Then we can define the discrete approximation for the exact eigenpair (λ,u) of
(3.1) based on the finite element space as:
Find (λ ,u ) ∈ R×V (Ω) such that b(u ,u ) = 1 and
h h 0h h h
a(u ,v ) = λ b(u ,v ), ∀v ∈ V (Ω). (3.2)
h h h h h h 0h
From (3.2), we know the following Rayleigh quotient expression for λ holds (see,
h
e.g., [2, 3, 23])
a(u ,u )
h h
λ = .
h
b(u ,u )
h h
Similarly, we know from [3, 7] the eigenvalue problem (3.2) has eigenvalues
0 < λ ≤ λ ≤ ··· ≤ λ ≤ ··· ≤ λ ,
1,h 2,h k,h Nh,h
6
and the corresponding eigenfunctions
u ,u ,··· ,u ,··· ,u ,
1,h 2,h k,h Nh,h
where b(u ,u ) = δ ,1 ≤ i,j ≤ N (N is the dimension of the finite element
i,h j,h ij h h
space V (Ω)).
0h
From the minimum-maximum principle (see, e.g., [2, 3]), the following upper
bound result holds
λ ≤ λ , i = 1,2,··· ,N .
i i,h h
Let M(λ ) denote the eigenspace corresponding to the eigenvalue λ which is
i i
defined by
M(λ ) = w ∈ V : w is an eigenvalue of (3.1) corresponding to λ
i i
(cid:8) and kwk = 1 , (3.3)
b
(cid:9)
where kwk = b(w,w). Then we define
b
p
δ (λ ) = sup inf kw−vk . (3.4)
h i 1
w∈M(λi)v∈V0h(Ω)
For the eigenpair approximations by the finite element method, there exist the
following error estimates.
Proposition 3.1. ([2, Lemma 3.7, (3.28b,3.29b)], [3, P. 699] and [7])
(i) For any eigenfunction approximation u of (3.2) (i = 1,2,··· ,N ), there is an
i,h h
eigenfunction u of (3.1) corresponding to λ such that ku k = 1 and
i i i b
ku −u k ≤ C δ (λ ).
i i,h 1,Ω i h i
Furthermore,
ku −u k ≤ C ρ (h)δ (λ ).
i i,h 0,Ω i Ω h i
(ii) For each eigenvalue, we have
λ ≤ λ ≤ λ +C δ2(λ ).
i i,h i i h i
Here and hereafter C is some constant depending on i but independent of the mesh
i
size h.
To analyze our method, we introduce the error expansion of eigenvalue by the
Rayleigh quotient formula which comes from [2, 3, 23].
Proposition 3.2. Assume (λ,u) is the true solution of the eigenvalue problem (3.1)
and 0 6= ψ ∈ H1(Ω). Let us define
0
a(ψ,ψ)
λ = .
b(ψ,ψ)
b
Then we have
a(u−ψ,u−ψ) b(u−ψ,u−ψ)
λ−λ = −λ .
b(ψ,ψ) b(ψ,ψ)
b
7
4 Multilevel local and Parallel algorithms
In this section, we present a new multilevel parallel algorithmto solve the eigenvalue
problem based on the combination of the local and parallel finite element technique
andthemultilevelcorrectionmethod. First, weintroduceanonecorrectionstepwith
the local and parallel finite element scheme and then present a parallel multilevel
method for the eigevalue problem.
For the description of the numerical scheme, we need to define some notation.
Given an coarsest triangulation T (Ω), we first divide the domain Ω into a number
H
of disjoint subdomains D , ···, D such that m D¯ = Ω¯, D ∩D = ∅ (see Figure
1 m j=1 j i j
2), then enlarge each Dj to obtain Ωj that aSligns with TH(Ω). We pick another
sequence of subdomains G ⊂⊂ D ⊂ Ω ⊂ Ω and (see Figure 2)
j j j
G = Ω\(∪m G¯ ).
m+1 j=1 j
In this paper we assume the domain decomposition satisfies the following property
Figure 2: m = 4
m
kvk2 . kvk2 (4.1)
ℓ,Ωj ℓ,Ω
Xj=1
for any v ∈ Hℓ(Ω) with ℓ = 0, 1.
4.1 One correction step
First, we present the one correction step to improve the accuracy of the given eigen-
value and eigenfunction approximation. This correction method contains solving an
auxiliary boundary value problem in the finer finite element space on each subdo-
main and an eigenvalue problem on the coarsest finite element space.
For simplicity of notation, we set (λ,u) = (λ ,u ) (i = 1,2,··· ,k,···) and
i i
(λ ,u ) = (λ ,u ) (i = 1,2,··· ,N ) to denote an eigenpair and its corresponding
h h i,h i,h h
8
approximationofproblems(3.1)and(3.2), respectively. Fortheclearunderstanding,
we only describe the algorithm for the simple eigenvalue case. The corresponding
algorithm for the multiple eigenvalue case can be given in the similar way as in [18].
In order to do the correction step, we build original coarsest finite element space
V (Ω) on the background mesh T (Ω). This coarsest finite element space V (Ω)
0H H 0H
will be used as the background space in our algorithm.
Assume we have obtained an eigenpair approximation (λ ,u ) ∈ R×V (Ω).
hk hk 0hk
The one correction step will improve the accuracy of the current eigenpair approxi-
mation (λ ,u ). Let V (Ω) ba a finer finite element space such that V (Ω) ⊂
hk hk 0hk+1 0hk
V (Ω). Here we assume the finite element spaces V (Ω) and V (Ω) are con-
0hk+1 0hk 0hk+1
sistent with the domain decomposition and V (Ω) ⊂ V (Ω). Based on this finer
0H 0hk
finite element space V (Ω), we define the following one correction step.
0hk+1
Algorithm 4.1. One Correction Step
We have a given eigenpair approximation (λ ,u ) ∈ R×V (Ω).
hk hk 0hk
1. Define the following auxiliary boundary value problem:
For each j = 1,2,··· ,m, find ej ∈ V (Ω ) such that
hk+1 0hk+1 j
a(ej ,v ) = λ b(u ,v )−a(u ,v ), ∀v ∈ V (Ω ). (4.2)
hk+1 hk+1 hk hk hk+1 hk hk+1 hk+1 0hk+1 j
Set uj = u +ej ∈ V (Ω ).
hk+1 hk hk+1 hk+1 j
2. Conestruct u ∈ V (Ω) such that u = uj in G (j = 1,··· ,m)
hk+1 0hk+1 hk+1 hk+1 j
and u = um+1 in G with um+1 being defined by solving the following
hk+1 e hk+1 m+1 hk+1 e e
problem:
e e e
Find um+1 ∈ V (G ) such that um+1| = uj (j = 1,··· ,m)
hk+1 hk+1 m+1 hk+1 ∂Gj∩∂Gm+1 hk+1
and
e a(um+1,v ) = λ b(u ,ve ), ∀v ∈ V e (G ). (4.3)
hk+1 hk+1 hk hk hk+1 hk+1 0hk+1 m+1
3. Define a neew finite element space V = V (Ω) + span{u } and solve
H,hk+1 0H hk+1
the following eigenvalue problem:
e
Find (λ ,u ) ∈ R×V such that b(u ,u ) = 1 and
hk+1 hk+1 H,hk+1 hk+1 hk+1
a(u ,v ) = λ b(u ,v ), ∀v ∈ V . (4.4)
hk+1 H,hk+1 hk+1 hk+1 H,hk+1 H,hk+1 H,hk+1
Summarize the above three steps into
(λ ,u ) = Correction(V (Ω),λ ,u ,V (Ω)),
hk+1 hk+1 0H hk hk 0hk+1
where λ and u are the given eigenvalue and eigenfunction approximation, re-
hk hk
spectively.
9
Theorem 4.1. Assume the currenteigenpairapproximation(λ ,u ) ∈ R×V (Ω)
hk hk 0hk
has the following error estimates
ku−u k . ε (λ), (4.5)
hk 1,Ω hk
ku−u k . ρ (H)ε (λ), (4.6)
hk 0,Ω Ω hk
|λ−λ | . ε2 (λ). (4.7)
hk hk
Then after one step correction, the resultant approximation (λ ,u ) ∈ R ×
hk+1 hk+1
V (Ω) has the following error estimates
0hk+1
ku−u k . ε (λ), (4.8)
hk+1 1,Ω hk+1
ku−u k . ρ (H)ε (λ), (4.9)
hk+1 0,Ω Ω hk+1
|λ−λ | . ε2 (λ), (4.10)
hk+1 hk+1
where ε (λ) := ρ (H)ε (λ)+ε2 (λ)+δ (λ).
hk+1 Ω hk hk hk+1
Proof. We focus on estimating ku−u k . First, we have
hk+1 1,Ω
ku−u k . ku−Pe uk +ku −P uk , (4.11)
hk+1 1,Ω hk+1 1,Ω hk+1 hk+1 1,Ω
and e e
m
ku −P uk2 = kuj −P uk2 +kum+1 −P uk2 . (4.12)
hk+1 hk+1 1,Ω hk+1 hk+1 1,Gj hk+1 hk+1 1,Gm+1
Xj=1
e e e
From problems (2.4), (3.1) and (4.2), the following equation holds
a(uj −P u,v) = b(λ u −λu,v), ∀v ∈ V (Ω ),
hk+1 hk+1 hk hk 0hk+1 j
for j = 1,2,··e· ,m. According to Proposition 2.3
kuj −P uk . kuj −P uk +kλ u −λuk
hk+1 hk+1 1,Gj hk+1 hk+1 0,Ωj hk hk −1,Ωj
. kuj −u k +ku −P uk +kλ u −λuk . (4.13)
ehk+1 hk 0,Ωj hk e hk+1 0,Ωj hk hk 0,Ωj
We will estimeate the first term, i.e. kej k by using the Aubin-Nitsche duality
hk+1 0,Ωj
argument.
Given any φ ∈ L2(Ω ), there exists wj ∈ H1(Ω ) such that
j 0 j
a(v,wj) = b(v,φ), ∀v ∈ H1(Ω ).
0 j
Let wj ∈ V (Ω ) and wj ∈ V (Ω ) satisfying
hk+1 0hk+1 j H 0H j
a(v ,wj ) = a(v ,wj), ∀v ∈ V (Ω ),
hk+1 hk+1 hk+1 hk+1 0hk+1 j
a(v ,wj ) = a(v ,wj), ∀v ∈ V (Ω ).
H H H H 0H j
10