Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic - - PowerPoint PPT Presentation
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems Svetozar Margenov Institute for Information and Communication Technologies, Bulgarian Academy of Sciences, Acad. G. Bonchev Str. Bl. 25-A, 1113 Sofia, Bulgaria Robust
Svetozar Margenov
Institute for Information and Communication Technologies, Bulgarian Academy of Sciences, Acad. G. Bonchev Str. Bl. 25-A, 1113 Sofia, Bulgaria Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 1/37
Collaboration: Ivan Georgiev, Johannes Kraus
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 2/37
Consider the weak formulation of a given elliptic b.v.p. in the form
a(u, v) = F(v), ∀v ∈ V,
and the related FEM problem
ah(uh, vh) = Fh(vh), ∀vh ∈ Vh.
Γ Γ
N D 1 2 3 4 5 7 8 9 10 11 12 13 14 15 16 17 18 19 20 6
We are interested in the efficient solution of the resulting large-scale FEM linear systems
Au = f.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 3/37
For large-scale problems, the iterative methods have advantages due to their better/optimal computational complexity and storage requirements. The Conjugate Gradient Method is the best iterative solution framework for large scale FEM systems. The construction of robust Preconditioned Conjugate Gradient (PCG) solution methods is addressed to some special properties of the stiffness matrix A, among which are that:
A is symmetric and positive definite (SPD); A is large and even very large but sparse.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 4/37
The number of iterations for the Conjugate Gradient (CG) method behaves as
n(CG)
it
= O
λmin .
Unfortunately, κ(A) = O(h−2), i.e., in 2D case κ(A) = O(N). To accelerate the convergence rate, the Preconditioned CG (PCG) method is used, for which
n(P CG)
it
= O
where the SPD matrix C is called preconditioner. Then, the efficient preconditioning strategy is determined by the conditions:
κ(C−1A) << κ(A) N(C−1v) << N(A−1v)
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 5/37
The preconditioner is called optimal, if
κ(C−1A) = O(1) N(C−1v) = O(N)
It is well known, that preconditioners based on various multilevel extensions of two-level FEM lead to iterative methods which have an optimal order computational complexity with respect to the size of the system. Unfortunately, the stiffness matrix A becomes additionally ill-conditioned when, e.g., the coefficients of the elliptic operator become more anisotropic or, equivalently, when the mesh aspect ratio increases. The condition number
κ(A) deteriorates also for various important parameter dependent problems
like almost incompressible elasticity, Navier-Stokes equations, etc. For such additionally ill-conditioned problems we need specially designed robust preconditioners.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 6/37
µ-FEM analysis of bone microstructure
The geometry of the solid phase of at a millimeter scale is obtained by a computer tomography (CT) image at a micron scale.
CT image of a trabecular bone specimen
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 7/37
µ-FEM analysis of coal-polyurethane geocomposite
The voxel data represent a coal-polyurethane material. The domain is cubic 75x75x75mm, and is non-uniform in all directions 35x110x110
Polyurethane [0.1, 0.25], E [200, 2100]MPa.
Figure: CT image of a coal-polyurethane geocomposite brick
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 8/37
Reservoir simulation The model has been originally generated for use in the PUNQ project (http://www.nitg.tno.nl/punq/index.htm). The fine scale model size is
1.122 × 106 cells.
Figure: Porosity for the whole model Figure: Lower fluvial part with clearly visible channels
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 9/37
Consider the SPD matrix A, written in the form:
A = A11 A12 A21 A22 = A11 A21 S I1 A−1
11 A12
I2 ,
where S states for the Schur complement
S = A22 − A21A−1
11 A12.
Let the additive (CA) and multiplicative (CM ) two-level preconditioners are defined as follows:
CA = A11 A22 , CM = A11 A21 A22 I A−1
11 A12
I .
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 10/37
The following estimates hold
κ
A A
1 − γ , κ
M A
1 1 − γ2
where γ ∈ [0, 1) stands for the constant in the strengthened Cauchy-Bunyakowski-Schwarz (CBS) inequality, corresponding to the block splitting of A. The presented condition number estimates clearly state the importance of such splittings where the CBS constant γ is uniformly bounded.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 11/37
Let us consider the sequence of nested triangulations T1 ⊂ T2 ⊂ . . . ⊂ Tℓ, corresponding FE spaces V1 ⊂ V2 ⊂ . . . Vℓ and related stiffness matrices
A(1), A(2), . . . , A(ℓ). The goal is to solve the finest discretization FEM system A(ℓ)u(ℓ) = f (ℓ).
Consider the 2 × 2 block presentation of A(k+1) corresponding to the splitting
A(k+1) = A(k+1)
11
A(k+1)
12
A(k+1)
21
A(k+1)
22
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 12/37
Let us consider the recursive multilevel generalization of the multiplicative two-level method introduced originally by Axelsson and Vassilevski (1989).
C(1) = A(1);
for k = 1, 2, . . . , ℓ − 1
C(k+1) = A(k+1)
11
A(k+1)
21
˜ A(k) I A(k+1)−1
11
A(k+1)
12
I ,
where the Schur complement approximation is stabilized by
˜ A(k)−1 =
A(k)−1.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 13/37
The acceleration polynomial is explicitly defined by
pβ(t) = 1 + Tβ 1 + α − 2t 1 − α
1 + α 1 − α
were α ∈ (0, 1) is a properly chosen parameter, and Tβ stands for the Chebyshev polynomial of degree β. Theorem: There exists α ∈ (0, 1), such that the AMLI preconditioner C = C(ℓ) has
complexity is O(N), if β satisfies the condition
4 > β > 1
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 14/37
For the model problem:
γ2 < 3/4 ⇒ β ∈ {2, 3} κ(A(k+1)
11
) = O(1) ⇒ NAMLI = O(N)
Construction of Robust Algebraic Multilevel Preconditioners: Uniform estimates of the CBS constant with respect to anisotropy and/or possible small parameters. Optimal order preconditioning (approximation) of the first pivot block
A(k+1)
11
with respect to anisotropy and/or possible small parameters.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 15/37
Consider the elliptic problem
−
2
∂ ∂xi
∂u ∂xj
in
Ω,
where Ω is a polyhedral domain, with proper boundary conditions on ∂Ω. Its variational formulation is : seek u ∈ H1
g(Ω) such that
pace2mm a(u, v) =
fv
for all v ∈ H1
0(Ω), where
a(u, v) =
2
aij ∂u ∂xi ∂v ∂xj ,
and where the function spaces H1
g(Ω) and H1 0(Ω) incorporate the Dirichlet
portion of the boundary conditions. Further, the matrix [aij] is assumed to be symmetric and positive definite (SPD).
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 16/37
Mesh anisotropy Coefficient anisotropy
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 17/37
The analysis for an arbitrary triangle (e) with coordinates (xi, yi), i = 1, 2, 3 can be done on the reference triangle (
e), with coordinates (0, 0) (1, 0) and (0, 1). Transforming the finite element function between these triangles, the
element bilinear form becomes
ae(u, v) = ae
e(
u, v) =
e
∂ u ∂ x, ∂ u ∂ y
(x2 − x1)(y2 − y1) (x3 − x1)(y3 − y1)
−1
a11 a12 a21 a22 (x2 − x1)(x3 − x1) (y2 − y1)(y3 − y1)
−1
∂ v ∂ x, ∂ v ∂ y T
where 0 <
x, y < 1, i.e. it takes the form ae
e(
u, v) =
e
∂ u ∂ xi ∂ v ∂ xj ,
and where the coefficients
aij depend on both triangle (e) and the coefficients aij in the differential operator.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 18/37
Corollary: It suffices for the analysis of uniform bounds for the finite element method to consider the reference triangle and arbitrary coefficients [aij], or alternatively, for the operator −∆ and an arbitrary triangle e. Given a coarse triangle, it can be subdivided in four congruent triangles by joining the mid-edge nodes, called the h−version. Alternatively, we can use piecewise quadratic basis functions in the added node points with support on the whole triangle, called the p−version.
(3) (5) (4) (1) (6) (2) Support for the piece-wise linear basis function ϕ
6
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 19/37
Theorem: For any finite element triangular mesh, where each element has been refined into congruent elements, it holds γ2
2 = 4
3γ2
1,
where γ1, γ2 are the CBS constants for the piecewise linear and piecewise quadratic finite elements, respectively. Theorem: Maitre and Musy [1982], Axelsson [1999], Blaheta, Margenov, Neytcheva [2004] The following estimate holds uniformly with respect to coefficient and mesh anisotropy
γ2
1 < 3
4
for both, conforming and nonconforming linear finite elements. Corollary:
γ2 < 1
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 20/37
Theorem: Axelsson, Padiy [1999], Axelsson, Margenov [2003] The additive preconditioner [A] of A11 has an optimal order convergence rate with a relative condition number uniformly bounded by
κ
−1A11
4(11 + √ 105) ≈ 5.31,
which holds independent on shape and size of each element and on the coefficients in the differential operator.
Connectivity pattern for the additive preconditioner [A].
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 21/37
Table: Number of PCG/GCG iterations for C.-R. FEM elliptic problem TP1 (a) TP1 (b) TP1 (c)
ℓ
V W2 W3 V W2 W3 V W2 W3 DA variant of the preconditioner 1 8 11 11 10 14 14 7 10 9 2 12 13 12 19 15 14 14 18 17 3 15 13 12 34 16 14 26 22 17 4 19 13 11 62 18 14 51 26 16 5 26 13 12 114 18 14 102 25 16 FR variant of the preconditioner 1 8 8 8 9 9 9 6 7 6 2 11 8 8 15 10 9 10 7 6 3 14 8 8 20 10 9 16 7 6 4 16 8 7 25 9 9 21 6 6 5 20 8 8 29 9 9 19 6 6 The robustness with respect to mesh anisotropy is tested for: (a) θ1 = 90o, θ2 = θ3 = 45o; (b)
θ1 = 156o, θ2 = θ3 = 12o; (c) θ1 = 177o, θ2 = 2o, θ3 = 1o.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 22/37
In the case of highly heterogeneous porous media the finite volume and mixed finite element methods have proven to be accurate and locally mass
enforced by Lagrange multipliers. Arnold and Brezzi have demonstrated that after the elimination of the unknowns representing the pressure and the velocity from the algebraic system the resulting Schur system for the Lagrange multipliers is equivalent to a discretization by Galerkin method using Crouzeix-Raviart linear finite elements. Further, such a relationship between the mixed and nonconforming finite element methods has been studied for various finite element spaces. Galerkin methods based on Crouzeix-Raviart and Rannacher-Turek finite elements has been also used in the construction of locking-free approximations for parameter-dependent problems.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 23/37
Rannacher-Turek nonconforming elements:
1 2 4 3 5 6
Figure 4. The reference element
For the variant MP (mid point), {ˆ
φi}6
i=1 are found by the interpolation condition
ˆ φi(bj
Γ) = δij, where bj Γ, j = 1, 6 are the centers of the faces of the cube ˆ
e: ˆ φ1 = (1 − 3x + 2x2 − y2 − z2)/6 ˆ φ2 = (1 + 3x + 2x2 − y2 − z2)/6 ˆ φ3 = (1 − x2 − 3y + 2y2 − z2)/6 ˆ φ4 = (1 − x2 + 3y + 2y2 − z2)/6 ˆ φ5 = (1 − x2 − y2 − 3z + 2z2)/6 ˆ φ6 = (1 − x2 − y2 + 3z + 2z2)/6
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 24/37
Alternatively, the variant MV (mean value) corresponds to the 3D integral mean-value interpolation condition |Γj
ˆ e|−1 Γj
ˆ e
ˆ φidΓj
ˆ e = δij , where Γj ˆ e are the
faces of the reference element ˆ
e: ˆ φ1 = (2 − 6x + 6x2 − 3y2 − 3z2)/12 ˆ φ2 = (2 + 6x + 6x2 − 3y2 − 3z2)/12 ˆ φ3 = (2 − 3x2 − 6y + 6y2 − 3z2)/12 ˆ φ4 = (2 − 3x2 + 6y + 6y2 − 3z2)/12 ˆ φ5 = (2 − 3x2 − 3y2 − 6z + 6z2)/12 ˆ φ6 = (2 − 3x2 − 3y2 + 6z + 6z2)/12
Figure 5. Uniformly refined macroelement
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 25/37
Kraus, Margenov, Robust Algebraic Multilevel Methods and Algorithms, de Gruyter [2009]
The multiplicative AMLI preconditioner is defined as follows,
H(0) = A(0), and for k = 1, . . . , ℓ, H(k) = B(k)
11
21
I(k)
1
B(k)−1
11
A(k)
12
I(k)
2
,
where
denotes that certain stabilization technique is performed. One particular stabilization is via a matrix polynomial, namely,
≡ S(k) = A(k−1) I − P k
νk(H(k−1)−1A(k−1))
−1 .
In this case, the following optimality condition for νk holds true:
1
where γ < 1 stands for the constant in the strengthened CBS inequality.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 26/37
Table 1: Multilevel behavior of γ2 for FR algorithm: 3D case
variant
ℓ ℓ − 1 ℓ − 2 ℓ − 3 ℓ − 4 ℓ − 5
MP 0.38095 0.39061 0.39211 0.39234 0.39237 0.39238 MV 0.5 0.4 0.39344 0.39253 0.39240 0.39238
MV MP 0.39238 0.38095 0.5 l l−5 l−4 l−3 l−2 l−1
Figure 6. Multilevel behavior of γ2 for FR algorithm: 3D case
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 27/37
Georgiev, Kraus, Margenov, Computing [2008]
The considered second order elliptic problem has a jumping coefficient a(e) which is elementwise constant.
Table: Linear AMLI W-cycle: number of PCG iterations: 3D tests MP
h−1
8 16 32 64 128 DA
ε = 1
9 10 10 10 10
ε = 10−3
9 10 10 10 10 FR
ε = 1
8 9 9 9 9
ε = 10−3
8 9 9 9 9 MV
h−1
8 16 32 64 128 DA
ε = 1
12 15 15 16 16
ε = 10−3
12 15 16 16 16 FR
ε = 1
10 12 12 12 12
ε = 10−3
10 12 12 12 12
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 28/37
For the case of “binary material”, the coefficient a(e) is initialized randomly, taking either of the values 1 or ε, where 1 occurs with a probability p. In the last case, the coefficient on each element is a uniformly distributed random number in (0,1).
Table 3: Non-linear AMLI W-cycle: GCG iterations for problem with random coefficients: 3D tests
p = 1/2
FR-MV
h−1
8 16 32 64 128
ε = 10−1
9 9 9 9 9
ε = 10−2
21 22 22 21 21
ε = 10−3
42 59 58 56 54
p = 1/10
FR-MV
h−1
8 16 32 64 128
ε = 10−1
9 9 9 9 9
ε = 10−2
17 22 22 22 22
ε = 10−3
29 60 55 50 50 Random coefficient: FR-MV
h−1
8 16 32 64 128
α(e) ∈ (0, 1)
22 39 43 34 35
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 29/37
The bone specimen is considered as a composite linear elastic material of solid and ”fluid” phases. The related coefficients are given by the bulk moduli
ks = 14 GPa and kf = 2.3 GPa, and the Poisson ratios νs = 0.325 and νf = 0.5. In our tests we vary νf ∈ {0.4, 0.45.0.49}.
Table 4: Convergence results for nonlinear AMLI W-cycle # voxels
νf = 0.4 νf = 0.45 νf = 0.49 163
14 19 46
323
14 20 50
643
14 21 53
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 30/37
Kraus, Margenov [2009], Georgiev, Kraus, Margenov [2010,2011], Kraus [2011]
Robust two-level results for the case of high-contrast and high-frequency (HC-HF) elliptic problems are presented in this section. The assumption is that the strong coefficient jumps can be resolved on the finest mesh only! Let us consider the symmetric interior penalty discontinuous Galerkin (IP-DG) finite element method: Find uh ∈ V such that
Ah(uh, v) = (f, v), ∀ v ∈ V,
where
Ah(uh, v) ≡ (a∇uh, ∇v)h + α
F [
[uh] ], [ [v] ]
− {a∇uh}, [ [v] ]F0∪FD −[ [uh] ], {a∇v}F0∪FD .
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 31/37
In the setting of IP-DG discretizations the construction of a face-based hierarchical splitting is proposed. The idea originates in the observation that the related global stiffness matrix can also be assembled from local matrices associated with the element faces.
Overlap of macro superelements and coarsening
Theorem: Let us consider a hierarchical basis two-level transformation based on interpolation with limiting weights. If α ≥ α0 = 25 then the estimate γ2
G ≤ 0.75
is uniform with respect to arbitrary jumps of the coefficients a(e) ∈ (0, 1]. Similar results were recently obtained for the case of bilinear conforming finite elements.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 32/37
Let us consider the IP-DG approximation in the form: Find uh ∈ Vh such that for all vh ∈ Vh the following relation holds:
A(uh, vh) = L(vh), A(uh, vh) = (a∇huh, ∇hvh)T −
[vh] ]
+θ
[uh] ], {a∇hvh}βe
F [
[uh] ], [ [vh] ]
where βe =
a− a++a− , ke := 2a+a− a++a− , and {au}βe = κe{u}.
The related IP-DG blinear form is continuous and coercive with constants independent of the mesh size h and coefficient a.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 33/37
The conditions to be met for optimal order multilevel methods, which are based
1 + γ 1 − γ < β < ̺, 1
Since the reduction factor of the number of DOF in our case is ̺ = 4, we can afford up two third degree stabilization (β = 3) in this setting. For the problem with jumps in the coefficient (a(x) = a(x)I), γ2 < 0.75. The multiplicative method will be of optimal order for β = 2. For the problem with anisotropic coefficient matrix γ2 < 0.8. The additive and the multiplicative methods will be of optimal order for β = 3.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 34/37
Table: Number of iterations and estimated condition number for the PCG: ǫ = 10−10, H=1/10, h=1/100, a = diag(η , 1), see the left figure; MS - multiscale coarse space, EMF - energy minimizing coarse space, and LSM - local spectral multiscale coarse space.
η
LIN MS EMF LSM (bilin. χi) LSM (MS χi)
103
113(1.48e+2) 122(1.51e+2) 115(1.81e+2) 53(23.21) 55(26.9)
104
257(1.35e+3) 258(1.28e+3) 231(9.70e+2) 41(53.63) 28(5.82)
105
435(1.34e+4) 483(1.26e+4) 416(9.64e+3) 28(5.642) 29(6.02)
106
627(1.34e+5) 709(1.27e+5) 599(9.63e+4) 30(5.753) 29(6.04) Dim 81=0.79% 81=0.79% 81=0.79% 732=7.19% 497=4.87% Efendiev, Galvis, Lazarov, Margenov, Ren [2011]
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 35/37
Discussions on coarse space dimension reduction: The strongly anisotropic channels cause a substantial increase of the size
To avoid this, we can replace the coarse solve RT
0 A−1 0 R0 by
RT A−1
0 R0 + RT anA−1 anRan.
The matrix
A0 is a small dimensional coarse matrix. The matrix Aan is
acting on the fine-mesh DOFs restricted to the high-anisotropy channels, locally constructed by preserving the strongest off-diagonal entries.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 36/37
TWO RELATED OPEN PROBLEMS: Robust preconditioning of higher order FEM anisotropic problems. High-contrast and high-frequency problems: from two-level to multilevel.
Robust Multilevel Methods for Anisotropic Heterogeneous Elliptic Problems – p. 37/37