Some numerical experiences on convergence criteria for iterative finite element solvers

Computers and Geotechnics 36 (2009) 1272–1284

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Some numerical experiences on convergence criteria for iterative finite element solvers X. Chen a,*, K.K. Phoon b a b

Department of Civil Engineering, National University of Singapore, Blk E2, #04-01, 5 Engineering Drive 2, Singapore 117576, Singapore Department of Civil Engineering, National University of Singapore, Singapore

a r t i c l e

i n f o

a b s t r a c t

Article history: Received 14 January 2009 Received in revised form 13 May 2009 Accepted 14 May 2009 Available online 12 June 2009 Keywords: Convergence criteria Iterative method Generalized Jacobi Modified SSOR Symmetric indefinite Decoupled relative residual norms

Several popular convergence criteria which are frequently used in practical finite element computations are investigated for two kinds of systems: the symmetric positive definite linear system and the symmetric indefinite system involving two distinct variables (displacement and pore fluid pressure). For the first system, the relative residual norm and the relative improvement norm are satisfactory as long as boundary fixities are handled appropriately. For the second system, the relative improvement norm must be adopted with greater care. It was further shown numerically that decoupled relative residual norms can be attractive alternates to the current global stopping criterion. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction In finite element analysis, we often need to solve a series of large-scale linear system of equations,

A 2 RNN

Ax ¼ b;

and x;

b 2 RN

N

ð1Þ NN

where R is the N-dimensional real Euclidean space and R is a real N  N matrix. To solve these large-scale linear equations, efficient linear solution techniques are needed because solving large-scale linear equations is often the most computationally expensive part in finite element analysis. During the past decades, two broad categories of linear solvers, namely direct solvers and iterative solvers, have been implemented in commercial softwares. The obvious difference between the two types of linear solvers is that iterative solvers only provide approximate solutions within a prescribed accuracy. Generally speaking, iterative methods are more suitable for large three-dimensional problems because for direct solution methods, the high memory requirement may limit their applications and the out-of-core facility may significantly slow down the computing speed. The crucial practical questions for iterative solvers are: when to terminate an iterative process and how to define an ‘‘acceptable” approximate solution. Stopping iterations pre-

* Corresponding author. Tel.: +65 6516 6783; fax: +65 6779 1635. E-mail addresses: [email protected] (X. Chen), [email protected] (K.K. Phoon). 0266-352X/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2009.05.012

maturely may lead to inaccurate solutions while prolonged iterations increase runtime without proportionate gain in accuracy. A robust convergence criterion that can deliver reliable solutions is of practical importance. Many convergence criteria have been proposed for linear equations arising from finite element computations. Some are based on the residual-norm, while others are based on the error-norm. From these two basic categories, some variants can be developed in terms of different vector norms such as Euclidean norm (i.e. 2norm), infinity norm (i.e. maximum norm) and energy norm (i.e. A-norm). Axelsson and Kaporin [1] gave some estimates on the lower and upper bounds for the PCG iteration error measured both in A-norm and 2-norm. The estimates can be used to assess the quality of the approximate solution of PCG or can be used as convergence criteria. Arioli et al. [2] proposed a residual dual-norm criterion for positive definite problems, but the difficulty is how to efficiently incorporate the evaluation of krk kH1 (here, H is the symmetric part of A, i.e. H = (A + AT)/2) in the iterative algorithms. Furthermore, they applied and extended their results to the indefinite saddle-point linear systems and nonlinear iterations arising from the mixed and mixed-hybrid finite element framework, but it is still limited to the global variable assessment [3]. Smith and Griffiths [4] used the relative ‘‘improvement” norm criterion for Krylov iterative methods. The obvious advantage for this criterion is that it only depends on the approximate solutions, which are outputs from any iterative algorithm. The relative residual norm criterion is another simple criterion, and it has been widely

1273

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

adopted for various applications because the current residual is usually the by-product during the iteration process. Recently, Picasso [5] proposed a convergence criterion based on a small fraction of the estimated error for nonlinear partial differential equations (PDEs) using anisotropic adaptive finite elements. Detailed descriptions of convergence criteria are given in Refs. [6,7]. It is quite accurate to say that most practitioners apply some popular criterion in conjunction with their solvers without quite knowing how accurate the ‘‘converged” solution is. Therefore, it is worth investigating the performance of some popular convergence criteria for practical applications in FEM codes. In this study, two types of symmetric linear systems of equations that are typical in engineering problems are examined. They are the symmetric positive definite (SPD) system and the symmetric indefinite system. The symmetric positive definite linear system is illustrated using a 2-D plane stress tensile plate analysis, while the symmetric indefinite system is illustrated using a 2-D plane strain consolidation analysis. The observations based on those 2-D problems are re-examined by using a larger 3-D footing problem with heterogeneous soils in drained and consolidation analysis, respectively. It has been recognized that convergence rate of a symmetric iterative method is mainly dominated by the spectrum of the preconditioned matrix, although the initial guess and the right-hand-side ‘‘force” vector can play some role. The effect of these influencing factors on the convergence behavior of the preconditioned symmetric Krylov iterative solver based on different convergence criteria are studied in this paper. It should be emphasized here that the present study only focuses on qualitative features of the convergence behavior that are related to the convergence criterion – not the actual rate of convergence such as iteration count. For the plate problem, we apply the displacement fixities using two different strategies: the penalty method and the variable elimination method. For the fully coupled consolidation analysis which generates a symmetric indefinite linear system, the difference between the convergence behaviors of symmetric iterative solver due to different convergence criteria is more apparent. In particular, when the preconditioned matrix has complex eigenvalues, it can produce local oscillations in the convergence curves that occasionally lead to premature termination for some criteria because of local minimum. In addition, it is also important to investigate if a single convergence criterion can be applied for different variable types, namely, the displacement variables and the excess pore pressure variables. The reason is that the solutions of some coupled problems contain different types of variables with significantly different magnitudes. The global based convergence criteria may be questionable and it may be necessary to examine the accuracy of displacements and excess pore pressures separately. The convergence of Krylov subspace solvers in finite precision arithmetic is generally more complex [8–10] than what is commonly portrayed in commercial softwares. Hence, numerical results presented in this paper are useful to practitioners. 2. Convergence criteria of iterative methods For a general boundary value problem (BVP) defined in domain X with Dirichlet boundary Cu and Neumann boundary Ct, the basic equations are given as 8 rrþb¼0 > > > > T 1 > > < e ¼ 2 ½ru þ ðruÞ  ð2Þ r ¼ De > > >  > u ¼ u on C u > > : t ¼ t on Ct

where r, e and u denote the stress field, the strain field and the displacement field, respectively. b is the body force vector, D is the  ; t are the fourth-order Cauchy tensor or constitutive tensor, and u prescribed displacement and the imposed traction, respectively. r is the gradient and r is the divergence operator. When applying finite element method to Eq. (2), we get the linear system of equation as

Ku ¼ f

ð3Þ

which is an expression of Eq. (1) in finite element framework. In Eq. (3), we have

(

R

P R BT DBdX ¼ e Xe BT DBdX P R P R f ¼ e Xe /ðxÞbdX þ e Ct /ðxÞðr  nÞdC K¼

X

ð4Þ

for the global assembled stiffness matrix and the external force vector, respectively. B = [/(x)] is the strain–displacement matrix, /(x) is the finite element basis function and n is the surface normal vector. 2.1. Krylov subspace iterative methods Krylov subspace iterative methods can be derived from different construction schemes of two essential ingredients: one is to construct the orthonormal basis vectors, and the other is to generate the approximate solutions. The Krylov subspace is constructed by a sequence of vectors as follows

Kk ðA; v Þ ¼ spanfv ; Av ; A2 v ; . . . ; Ak1 v g;

k ¼ 1; 2; . . .

ð5Þ

in which, Kk (A, v) is called the kth Krylov subspace of R generated by A w.r.t v, and usually v is taken as the initial residual r0. By making use of the sequence of vectors defined in Kk (A, v), we can construct a new sequence of basis vectors with orthonormal property. For symmetric matrix, the construction algorithm is called Lanczos method, while for nonsymmetric matrix the construction algorithm is called Arnold method. The approximate solution is constructed based on the initial guess solution x0 and the current Krylov subspace, Kk (A, v), in other words, N

xk 2 x0 þ Kk ðA; v Þ;

k ¼ 1; 2; . . .

ð6Þ

The approximation of xk is constructed based on different criteria such as the minimization criterion (e.g. [11–14]). Because the appropriate Krylov iterative solver depends on the problem, it is of practical interest to investigate its convergence based on different criteria for different problems. In particular, for nonsymmetric problems, different iterative solver exhibits quite different convergence behavior based on the same convergence criterion [15,16]. Although the present study only focuses on symmetric linear systems of equations, it does highlight that convergence criteria are not straightforward, in contrary to popular belief. Previous studies [7] showed that the same convergence criterion produces similar behavior in different symmetric solvers if the zero initial guess is given. For the SPD linear systems of equations, PCG method is apparently the most robust one, and thus it is employed in this study. Unlike the SPD linear systems, some alternative options such as MINRES, SYMMLQ and SQMR can be used for the symmetric indefinite linear systems of equations (e.g. [17]). However, among these methods, only SQMR solver allows for an arbitrary symmetric preconditioner (e.g. [18]). Usually, a suitable preconditioner should be adopted to accelerate the convergence of an iterative method. Preconditioning is equivalent to transforming the original linear system as 1 1 ~~ ~ P1 L AP R P R x ¼ P L b or Ax ¼ b

ð7Þ

~¼ where P = PLPR is the applied preconditioner and A is the ~ can be preconditioned matrix. The eigenvalue distribution of A 1 P1 L AP R

1274

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

regarded as an important indicator for the convergence rate. In the ~ is an identity extreme, if eigenvalues are clustered fully at 1, then A matrix and solution is determined in one step. Refer to [19,20] for more preconditioning methods of iterative solutions. 2.2. Convergence criteria In practice, to attain an acceptable solution, we have to terminate an iterative method using convergence criteria. In practical applications, two or more criteria can be combined in case one of them is ineffective. The convergence criteria for an iterative method in finite element analysis can be broadly categorized into the error-norm based criteria and the residual-norm based criteria. Different vector norms can be taken into account. The relative error norm criterion is:

R¼

kek k kx  xk k ¼  tol; ke0 k kx  x0 k

k ¼ 1; 2; . . . ; maxit

ð8Þ

pffiffiffiffiffiffiffiffiffi in which the 2-norm is often adopted, that is, kv k2 ¼ v T v for any vector v. maxit and tol are user specified maximum iteration number and the convergence tolerance, respectively. Note that R requires the exact solution (x) to be known apriori. Although it is not applicable, it serves as the ‘‘ideal” criterion to benchmark other more practical criteria in this study. In other words, the ‘‘best” criterion is one that tracks R as closely as possible. We introduce several popular convergence criteria in finite element analysis below: (1) Relative residual convergence criterion. The relative residual convergence criterion has been widely adopted (e.g. [6,11,12,21]), and it is defined as

Rr ¼

kr k k2 kb  Axk k2 ¼  tol; kr0 k2 kb  Ax0 k2

k ¼ 1; 2; . . . ; maxit

kA1 k2  kr 0 k2 kA1 k2  kr 0 k2 Rr with 1 ke0 k2 ke0 k2

ð10Þ

ð11Þ

From Eq. (11), we know that the Rr convergence curve may lie above or below the R curve, which will be demonstrated by two numerical examples below. In other words, Rr does not systematically underpredicts or overpredicts R. For a SPD matrix it is well known that [1]

kx  xk k2 krk k2  jðAÞ kxk2 kbk2

ð12Þ

where j() denotes the spectral condition number of a matrix, and it is computed by

jðAÞ ¼ kAk2  kA1 k2 ¼ kmax =kmin

kxk  xk1 k1  tol; kxk k1

k ¼ 1; 2; . . . ; maxit

ð14Þ

in which ||||1 represents the infinity norm. The relative ‘‘improvement” norm criterion has been applied and studied by some researchers (e.g., [4]). In fact, when the initial solution is set as x0 = 0, the relative ‘‘improvement” norm criterion is closely related to the relative residual convergence criterion, but this relation is also dependent on the property of matrix A. Theoretical analysis of the relation between Rr and Ri is difficult. Hence, numerical evaluation becomes important.

3. Problem associated with prescribed boundary fixities 3.1. Neumann boundary condition In Eq. (4), we only imposed the Neumann boundary conditions on Ct. In a finite element analysis, the Neumann boundary condition is also called as force boundary condition which can be imposed as follows

r  n ¼ t or

ð9Þ

We can rewrite Eq. (10) to make explicit the relation between R and Rr

R

Ri ¼

rij nj ¼ ti

ð15Þ

in which r = {r11, r22, r33, r12, r23, r31}T for 3-D problems and n is the surface normal vector with the form

Note that the residual norm, ||b  Axk||2 = ||Ax  Axk||2. While we cannot calculate the error norm ||x  xk||2 without knowledge of the exact solution (x), it is possible to calculate ||Ax  Axk||2 without x. Hence, there is an advantage to use the residual, although it is an indirect measure of the error. Based on the relation ek = A1rk, we know that the criterion has the following error bound,

kek k2  kA1 k2  krk k2 6 tol  kA1 k2  kr0 k2

pared numerically in Section 5. This potential difficulty in comparing criteria based on different indirect measure of ‘‘error” and different norm applies to other criteria discussed below. (2) Relative ‘‘improvement” norm criterion. The relative ‘‘improvement” norm criterion is defined based on the approximate solution as

ð13Þ

Here kmax and kmin are the largest and the smallest eigenvalues of A, respectively. If the initial guess is set as zero, Eq. (12) shows that the relative error norm is bounded by the relative residual norm multiplied with j(A), i.e. R 6 j(A)Rr. It is important to note that R and Rr can differ by orders of magnitude for very large j(A) (ill-conditioned matrices) when these criteria are being com-

2

n1

6 n¼4 0

0

0

n2

0

n2

0

n1

n3

0

0

n3

0

n2

n3

3

7 05 n1

ð16Þ

with ni (i = 1, 2, 3) are cosines of the boundary surface normal with respect to the i axis. 3.2. Dirichlet boundary condition The Dirichlet boundary condition is also called the displacement boundary condition in solid mechanics. Usually, we have two approaches to apply the prescribed displacement boundary condition. One approach is the penalty or large number method, and the other one is the variable elimination approach. Though the two approaches are widely known, they are still briefly reviewed to support subsequent discussions related to some convergence criteria. 3.2.1. Penalty/large number method Penalty method is widely used to assign the prescribed boundary conditions in many finite element packages such as the one in Ref. [22]. When applying a prescribed displacement ui = ci to the linear equation (3), we have

2

K 11 6K 6 21 6 . 6 . 6 . 6 6K 6 i1 6 . 6 . 4 .

K 12



K 1i



K 22 .. .

 .. .

K 2i .. .

 .. .

K i2 .. .

 .. .

K ii .. .

 .. .

K N1

K N2

   K Ni

3 2 3 f1 u1 7 6 6 7 K 2N 76 u2 7 7 6 f2 7 7 6 7 6 .. 7 76 . 7 6 . 7 . 76 .. 7 6 .. 7 76 7¼6 7 7 6 7 6 K iN 7 76 ui ¼ ci 7 6 fi 7 7 6 7 6 .. 7 76 .. 7 6 .. 7 . 54 . 5 4 . 5 K 1N

   K NN

32

uN

fN

ð17Þ

1275

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

The penalty method can apply the boundary condition, but retains the number of rows (or columns) of the stiffness matrix. In this method, a large number such as w = 1E+40 is used to modify the stiffness matrix as well as the right hand side (RHS) vector in the following steps, (1) Replace the diagonal entry Kii with w; (2) Replace the corresponding component fi in RHS with w  ci; (3) However, the diagonal large number may artificially introduce the ill-condition to the stiffness matrix. One possible remedy is to scale the ith row and the ith column, respectively, before solving this linear algebraic equation. One pffiffiffiffi can scale the ith row and ith column by 1= w so that the diagonal entry becomes unit. Correspondingly, the ith compffiffiffiffi ponent in RHS should also be scaled to be w  ci . It should be noted that once the linear algebraic equation is solved, pffiffiffiffi the solution should be scaled back by 1= w. When Eq. (3) is restrained with the penalty method as described above, the resultant linear equation is denoted as

K1 u1 ¼ f 1

ð18Þ

If the relative residual criterion is employed, however, there are still pffiffiffiffi some large numbers w  ci in the RHS. Direct application of this criterion may fail so that these components have to be removed, which leads to the following natural modification,

^ r ¼ k^r k k2  tol; R k^r 0 k2

k ¼ 1; 2; . . . ; maxit

k ¼ 1; 2; . . . ; maxit

(1) size(K2) 6 size(K1) = size(K); (2) In Eq. (21), the RHS vector contains some large components.

4. Coupled multi-phase problems Some typical multi-phase problems may involve different types of variables in their solutions. For instance, the fully coupled soil consolidation problem in geomechanics formulated by Biot [23] and the incompressible Stokes flow system [24] are such cases. In the linear system of equations arising from Biot’s consolidation equations, the displacement unknowns are coupled with the excess pore pressure unknowns. The coupled incremental Biot’s formulation is given as (e.g. [4,25–27])



ð19Þ

^ ¼ v :T for any vector v and T denotes the set of removed in which v variables. Similarly, when the fixities are applied by the penalty method the relative ‘‘improvement” norm criterion should also be modified, i.e.,

^ i ¼ k^xk  ^xk1 k1  tol; R k^xk k1

 ¼ uT and T is the variable set with where, K2 = K:-T,:T, u2 = u:T, u prescribed values. When assembling each new element stiffness matrix, the rows and columns of the element stiffness matrix are assembled into K2 if u R T; otherwise element-level elimination should be conducted in terms of RHS of Eq. (23). Therefore, the storage for the assembled matrix can be reduced at the cost of elimination steps. Obviously, by using the elimination method the dimension of the stiffness matrix is reduced. The main difference between Eqs. (18) and (23) can be summarized as

ð20Þ

K LT

K j1 u1 þ    þ K j;i1 ui1 þ K j;iþ1 uiþ1 þ    þ K jN uN ¼ fj  K ji ui ;

ðj ¼ 1; . . . ; N; j–iÞ

ð21Þ



 ¼

Df nþ1 DtQpex n

ð22Þ

Here :(i) denotes the variable set excluding the ith variable. To eliminate all equations with the prescribed values, we obtain

ð23Þ



ðn ¼ 1; 2; . . .Þ

ð24Þ

To simplify the notations, Eq. (24) is rewritten as



K L

T

L C



xu x



p

¼

  f g

ð25Þ

For the coupled stiffness matrix in Eq. (25), it has m (dimension of matrix K) positive eigenvalues and n (dimension of matrix C) negative eigenvalues [30]. Dependent on the preconditioner type, the preconditioned matrix may have purely real eigenvalues or complex eigenvalues. Applying the block Gaussian elimination to Eq. (25) results in the decoupled linear systems



Kxu ¼ f  Lxp Sxp ¼ LT K 1 f  g

K:ðiÞ;1 u1 þ    þ K:ðiÞ;i1 ui1 þ K:ðiÞ;iþ1 uiþ1 þ    þ K:ðiÞ;N uN

 K2 u2 ¼ f 2 ¼ f  K:T;T u

Dunþ1 Dpex nþ1

(1) Case b1: Df nþ1 –0; pex n ¼ 0 for load driven consolidation process, which happens only for the first load step; (2) Case b2: Df nþ1 ¼ 0; pex n –0 for time driven consolidation process, which happens when no external loads are applied; (3) Case b3: Df nþ1 –0; pex n –0 for concurrent load and time driven consolidation process, which may happen when additional loads are continuously added during the consolidation process.

In matrix form, Eq. (21) can be rewritten as

¼ f :ðiÞ  K:ðiÞ;i ui



where K, Q are the effective solid and fluid stiffness matrices, respectively, and C = hDtQ. Dt is the time increment, while h is the time interpolation parameter with h = 1 for the fully implicit scheme. Dun+1 and Dpex nþ1 are the incremental displacement and the incremental excess pore fluid pressure, respectively, at the nth time step. Dfn+1 is the incremental external forces. For Dt > 0, three right hand side (RHS) vectors can arise:

In the penalty approach, a large number is used and usually the resultant system is ill-conditioned. An alternative method which can obviate this problem is to set the diagonal entry as one and the off-diagonal entries as zeros, and replace the corresponding component fi in RHS with ci. Compared to the penalty approach without scaling, this method involves some initialization work. 3.2.2. Elimination method The penalty method is convenient to implement in a computer code. Any row or column index of the global stiffness matrix can be located quickly. By using the penalty method, however, we are solving some unknowns which we actually know. Therefore, the elimination method is proposed to obviate the unnecessary work and reduce the storage requirement. For example, the variable elimination method is employed in the finite element codes of [4]. The idea is that once we know the prescribed variable value xi = ci for the ith row equation, we can eliminate the corresponding ith columns in the other rows, which can be formulated as

L C

ð26Þ

in which S = C + LTK1L is the Schur complement matrix, and the corresponding residuals at the kth iteration can be given as

(

suk ¼ ðf  Lxpk Þ  Kxuk spk ¼ LT K 1 ðf  Lxpk Þ  Cxpk  g

ð27Þ

1276

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

in which we use the notations of suk and spk for the decoupled residuals to distinguish from the separated residuals, ruk and r pk , which are directly extracted from the global residual rk as



ru r k ¼ kp rk

"

 ¼

f  ðKxuk þ Lxpk Þ

#

ð28Þ

g  ðLT xuk  Cxpk Þ

In short, the decoupled residuals are determined by separating the global linear system to two systems, while the separated residuals are determined by separating the global solution vector into two vectors. Mathematically, suk ¼ ruk but spk –rpk . When the relative residual norm criterion is applied to the decoupled residuals, respectively, we have

8 ksu k u > < Rs ¼ ksku k22  tolu 0 p

> : Rps ¼ kspk k2  tolp ks k

with k ¼ 1; 2; . . . ; maxit

ð29Þ

0 2

in which suk and spk are calculated based on Eq. (27). If the relative residual norm criterion is constructed by using the separated residuals, r uk and r pk , as defined by Eq. (28), the resultant criteria are denoted as Rur and Rpr . Note that mathematically Rur is equal to Rus as highlighted above. Similarly, the separated relative ‘‘improvement” norm criteria are defined as:

8 kxu xu k1 u >  tolu < Ri ¼ kkxuk1 k1 p

0 p

with k ¼ 1; 2; . . . ; maxit

> : Rpi ¼ kxk xpk1 k1  tolp kx k

ð30Þ

0 1

Note that the separated/decoupled residual vectors are shorter than the global vector. In particular, for typical geotechnical engineering problems, the excess pore water pressure degrees of freedom constitute about 10% of the total degrees of freedom. 5. Numerical evaluation of convergence criteria Two types of linear systems of equations that are typical in engineering problems are examined using three numerical examples. They are the symmetric positive definite system and the symmetric indefinite system. First, we examine the convergence criteria for the symmetric positive definite linear system stemming from the tensile plate in plane stress analysis. The finite element mesh is composed of 8-node linear strain quadrilateral elements

with uniform size. The geometry dimension, material parameters of the plate and the prescribed displacement are given in Fig. 1. In this figure, E0 is the effective Young’s modulus, m is the Poisson ratio. In the second example, we study the convergence behavior of an iterative solver for a plane strain pile-group consolidation problem constructed from 8-node solid elements coupled with 4node fluid elements. Consolidation around a pile group has important negative skin friction ramifications. The model material zones and geometry details with the consideration of symmetry are generated by GeoFEA [29] as shown in Fig. 2. Apart from the parameters E0 and m for the solid effective stiffness of two soil layers and the pile material, we also need the hydraulic conductivities kx and ky in the x and y direction, respectively, for the fluid stiffness. For the symmetric indefinite formulation arising from fully coupled consolidation problems, SQMR solver has shown good performance and thus has been employed in many studies [26–33]. To accelerate the convergence of SQMR, we may use the recently developed modified SSOR (MSSOR) preconditioner which has the following preconditioning format [32],

PL ¼

LA þ

^ D

x

! ;

PR ¼

^ D

x

!1 UA þ

^ D

!

x

in which LA and UA are the strictly lower and upper triangular parts ^ is constructed by using the Generalized of matrix A, respectively, D Jacobi matrix [33] and x is the relaxation parameter. We use x = 1.0, which corresponds to the symmetric Gauss–Seidel (SGS) version, in the following studies. The third example is a 3-D footing problem as shown in Fig. 3, for which drained and consolidation analysis are carried out, respectively. The resultant symmetric positive definite and symmetric indefinite linear systems are used to further examine the observations from the above 2-D examples. Numerical experiments show that convergence behaviors (at least for symmetric matrices) are closely related to the eigenvalue distribution/clustering of the preconditioned matrix [34,35]. Our empirical observations are summarized as: (1) The convergence rates are dominated by the eigenvalue clustering; (2) Dramatic local oscillations of the convergence curves appear to be introduced by the imaginary part of the eigenvalues.

δy=0.001m

Γu2 Ω Γt1

Γt2

ly = 2m

E’ = 20000 MPa ν = 0.2 lx = 1m Γu1 Fig. 1. Plate subjected to the prescribed displacement.

ð31Þ

Fig. 2. Pile group in plane strain analysis.

1277

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

10

1

10

-2

R Ri

10

-5

Rr

10

-8

||rk||2

10

-11

10

-14

10-17 10-20

(a) 0 10

Fig. 3. Footing foundation modeled with two different soil layers.

1000

2000

1

10

-2

R Ri

10

-5

Rr ||rk||2

10-8 The above observations are illustrated by the examples below. It will be shown that observation (2) can pose a serious practical problem to proper application of convergence criteria.

10

-11

10

-14

10

-17

5.1. Plane stress plate example

10

-20

The purpose of this example as shown in Fig. 1 is to illustrate the problems associated with prescribed boundary fixities. 5.1.1. Convergence behaviors based on different criteria The elimination approach is first adopted to impose boundary fixities. The methods considered are Jacobi preconditioned CG (JCG), Jacobi preconditioned SQMR (JSQMR), and the 0-level incomplete Cholesky preconditioned CG [ICCG(0)] [36,37]. Note that for a symmetric positive definite preconditioner, SQMR is mathematically equivalent to MINRES so that the residual norms are truly minimized [15]. The ideal criterion, R in Eq. (8), can not be evaluated directly since the exact solution is unavailable. We evaluate the exact solution x separately using a direct solver. The performances of other criteria are measured against R. We have remarked in Section 2.2 the potential difficulty in comparing criteria based on different indirect measure of ‘‘error” and different norm. In Fig. 4a, it can be seen that Ri and Rr exhibit very similar convergence behaviors for the preconditioned CG methods and they do track R quite closely. For SQMR, the overall convergence trends for Ri and Rr are similar as well though local departures may occur (Fig. 4b). For Fig. 4 and other similar figures below, a horizontal line indicating the customary prescribed tolerance of 1E6 is plotted. The iteration counts corresponding to Ri and Rr are shown by the downward arrows. The values of R corresponding to these iteration counts are shown by the leftward arrows. It is noteworthy that the apparent residuals Ri and Rr are smaller than the true relative error R. Hence, for a prescribed tolerance, Ri and Rr will cause the solver to terminate before R. The true residual norm, i.e. ||rk||2, is also provided. Rr can be regarded as a scaled true residual norm (Eq. (9)). Because the initial guess may influence the convergence behavior of an iterative solver [38], Fig. 5 investigates three different cases of initial guess x0, namely: x0 = zeros(N, 1), x0 = ones(N, 1) and x0 = rand(N, 1) (MATLAB conventions for a vector with zeros, ones, and random numbers between 0 and 1, respectively). It can be seen that the initial guess may significantly influence the convergence behaviors of an iterative solver. For non-zero initial guesses, the

3000

(b) 0

1000

2000

101 10

R Ri

-2

Rr

10-5 10

3000

||rk||2

-8

10

-11

10

-14

10-17 10

-20

(c) 0

500

1000

1500

Fig. 4. Convergence behaviors of preconditioned iterative solvers based on different criteria for the plane stress plate problem: (a) JCG; (b) JSQMR; (c) ICCG(0).

convergence curves based on Rr now lie below the one based on Ri so that Ri should be a more conservative criterion. In Fig. 6, the eigenvalue distributions of preconditioned matrices using the 20  20 finite element mesh for the tensile plate are plotted. The eigenvalues of preconditioned matrices are distributed on the real axis (i.e. zero imaginary part), but the IC(0) preconditioned matrix (Fig. 5a) possesses a more clustered real eigenvalue distribution. The condition number [max(ki)/min(ki)] for the Jacobi preconditioned matrix and the IC(0) preconditioned matrix are 2.13E+05 and 3.80E+04, respectively. The convergence rates for symmetric matrix systems appear to be dominated by the eigenvalues of the preconditioned systems, and the spectral condition number is an effective measure of the eigenvalue distribution [34,35]. According to Benzi [19], good spectral properties which may lead to rapid convergence are: eigenvalues clustered around 1, less clustered points and small clustered radius. Referring to Figs. 4 and 6, it can be seen that the eigenvalues of Jacobi

1278

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

104 101 10-2

Rr

10

-5

10

-8

10

Imag(eig)

1

R Ri ||rk||2

max(λi) = 3.7E+00 min(λi) =1.74E−05

-1

(a) 0

1

2

3

4

-11

1

Imag(eig)

10-14 10-17 10

0

-20

(a) 0

500

1000

max(λi) = 2.52E+00 min(λi) = 6.64E−05

-1

(b) 0

104

R Ri

101 10-2

Rr

10-5 10

1500

0

1

2 Real(eig)

3

4

Fig. 6. Eigenvalue distribution of preconditioned matrices using the 20  20 finite element mesh for the plane stress plate problem: (a) Jacobi preconditioned matrix; (b) IC(0) (0-level Incomplete Cholesky) preconditioned matrix.

||rk||2

-8

convergence behavior – a slow linear convergence followed by a fast superlinear convergence before the plateau stage – can be explained from the perspective of eigenvalue distribution. The slow linear convergence is associated with isolated and extreme eigenvalues. The fast superlinear convergence is associated with clustered eigenvalues.

10-11 10-14 10-17 10-20

(b) 0

500

1000

1500

104 R Ri

101 10-2

Rr

10-5 10

||rk||2

-8

10-11 10-14 10-17 10-20

(c) 0

500

1000

1500

Fig. 5. Convergence behaviors of ICCG(0) based on global criteria for the plane stress plate problem: (a) x0 = zeros(N, 1); (b) x0 = ones(N, 1); (c) x0 = rand(N, 1).

preconditioned matrix and the IC(0) preconditioned matrix are all real and there are no dramatic local oscillations in their convergence curves. This provides some numerical evidence for observation (2) in Section 5. From Figs. 4 and 5, we further note that R exhibits a final plateau. The residual norm ||rk||2 will follow the same plateau if the residual, rk = b  Axk, is computed directly. However, it is not efficient to calculate rk directly because of the matrix–vector product Axk. In the PCG algorithm, the residual is updated more efficiently using rk = rk1  ak1qk1 because qk1 has already been computed. The directly computed rk and the updated rk are almost the same until the final plateau is reached. Thereafter, the updated rk is probably increasingly less accurate. However, the practical tolerance is usually set at a level that is much higher than the final plateau. Hence, the more efficient updated rk is acceptable in practice. Figs. 4 and 5 support the two-phase convergence behavior observed by some researchers [5,39]. Theoretically, the two-phase

5.1.2. Convergence problems associated with the penalty approach Table 1 compares the iteration count and true error as measured by R for the penalty approach and the elimination approach. The solver is JCG and iterations are terminated at a tolerance of 1E6 based on different convergence criteria. In this study, the spatial domain is discretized into a series of uniform mesh sizes from 5  5  5 to 50  50  50. N denotes the total number of equations (or DOFs) and iters denotes the required iterations to ‘‘converge”. For the elimination approach, the modified criteria are the same as the standard ones. For the penalty approach, the differences are pronounced as expected. When the standard criteria are applied, both Ri and Rr caused extremely early termination as highlighted by the asterisks. When measured with respect to R, we know that the solutions are not sufficiently accurate. For Rr, this premature termination is caused by some large numbers in the RHS which make the denominator ||r0||2 large. For Ri, this premature termination is caused by ||xk||2, which is very large in the first few iterations. In fact, these constrained equations corresponding to the fixed degrees of freedoms (DOFs) are decoupled from the rest of

Table 1 Convergence behaviors of JCG method based on different criteria with the tolerance 1E6. Mesh

N

Rr

^r R

Elimination method, iters/R 169 109/9.7E7 52 102 639 231/1.7E6 202 2479 464/3.2E6 502 15,199 1154/4.3E6

Same Same Same Same

Penalty method, iters/R 192 1a/4.8E3 52 682 1a/9.9E3 102 2562 1a/2.0E2 202 15,402 1a/5.0E2 502

125/1.2E8 265/1.4E8 539/3.4E8 1349/4.6E8

a

Denotes premature termination.

as as as as

Rr Rr Rr Rr

Ri

^i R

104/3.7E6 227/2.7E6 457/4.5E6 1133/1.0E5

Same Same Same Same

2a/4.6E3 2a/9.6E3 2a/2.0E2 2a/5.0E2

110/4.7E9 228/2.7E8 458/9.1E8 1135/4.9E7

as as as as

Ri Ri Ri Ri

1279

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

the DOFs. Hence, it is more reasonable to neglect these constrained ^ i are defined in Eqs. (19) ^ r and R equations. The modified criteria R ^ i indeed can produce accurate solutions and (20). From Table 1, R at iteration counts comparable to those from the elimination ap^ r , the iteration counts are higher, ranging from about proach. For R 14% higher for the 5  5  5 mesh to 19% higher for the 50  50  50 mesh. In contrast, for the elimination approach, Rr only produces iteration counts that are less than 5% higher than ^ i and R ^ r are those produced by Ri. Overall, the modified criteria R successful in preventing premature termination of the iterative solver and should be adopted when the penalty approach is used to impose boundary fixities. It should be mentioned that the residual norm criterion ||Axk  b||2 can effectively obviate the above premature termina-

tion associated with the penalty approach. However, convergences of iterative solvers based on the residual norm criterion ||Axk  b||2 are RHS dependent. For example, small RHS norms may lead to premature terminations of iterative solvers if the residual norm criterion is employed. Hence, the residual norm criterion is not appropriate for general problems. 5.2. Plane strain pile-group example The purpose of this example as shown in Fig. 2 is to demonstrate the reliability of different convergence criteria and to illustrate the effect of different types of variables in the solution vector.

10

4

10

1

R Ri

10

-2

Rr

10

-5

10-5

10

-8

10-8

10-11

10-11

10

104

-14

(a) 0

10

8000

12000

104

16000

R Ri

101

-14

(a) 0 10

4

10

1

Rr 10

-2

10

-2

10

-5

10

-5

10-8

10

-8

10-11

10-11

10-14

10-14

(b) 0 10

4

10

1

400

800

1200

1600 R Ri Rr

(b) 0

400

800

1200

1600

10

-2

Rr

-5

10

-8

10

-8

10

-11

10

-11

10

-14

10

-14

Fig. 7. Convergence behaviors of SQMR preconditioned by different preconditioners based on global criteria for the plane strain pile group problem (b1 RHS case): (a) GJ; (b) Pc; (c) MSSOR.

Rr

R Ri

10

3200

R Ri

1

-5

2400

12000

10

10

1600

8000

4

-2

800

4000

10

10

(c) 0

Rr

10-2

10

4000

R Ri

1

(c) 0

800 1600 2400 3200 4000 4800

Fig. 8. Convergence behaviors of SQMR preconditioned by different preconditioners based on global criteria for the plane strain pile group problem (b2 RHS case): (a) GJ; (b) Pc; (c) MSSOR.

1280

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

104 1

R Ri

10-2

Rr

10

10

-5

10

-8

10

-11

10

-14

(a) 0 10

4

10

1

4000

8000

12000

R Ri

5.2.1. Global convergence behaviors based on different criteria The solution vector contains two types of variables: displacement and the excess pore fluid pressure. Using a zero initial guess, Figs. 7–9 compare the convergence behaviors of SQMR preconditioned by three preconditioners, namely GJ, Pc (block constrained preconditioner in [40]) and MSSOR [32] for three RHS cases mentioned in Section 4, respectively. In contrast to the plane stress tensile plate example (Fig. 4), the convergence behaviors of SQMR based on Ri and Rr are not identical. It is evident that Ri is tracking R from below while Rr is tracking R from above before the final plateau of R. Hence, when the same tolerance is prescribed, the approximate solution obtained by using Ri can be significantly less accurate than that obtained by using Rr. More extensive results (iterations and relative error at the converged point) for three mesh sizes are reported in Table 2. The results clearly show that Ri criteria may produce significantly larger errors than the relative residual criteria for the specified tolerance 1E6. The relative errors due to the relative ‘‘improvement” norm criteria are usually

Rr

10

-5

10

-8

Imag(eig)

10-2

10-11

(a) 0

-14

(b) 0 10

400

800

1200

1600

4

R Ri

101 -2

10

-5

10

-8

10

-11

10

-14

(c) 0

Rr

1600

2400

1

2

3

4

5

max(λi) = 3.809E+00 min(λi) = 3.885E−06 1

2

0.3 0.2 0.1 0 -0.1 -0.2 -0.3

(c) 0 800

max(λi) = 4.147E+00 min(λi) = 3.885E−06

0.3 0.2 0.1 0 -0.1 -0.2 -0.3

(b) 0 Imag(eig)

10

Imag(eig)

10

0.3 0.2 0.1 0 -0.1 -0.2 -0.3

3

4

5

max(λi) = 1.000E+00 min(λi) = 3.893E−06 1

2 3 Real(eig)

4

5

3200

Fig. 9. Convergence behaviors of SQMR preconditioned by different preconditioners based on global criteria for the plane strain pile group problem (b3 RHS case): (a) GJ; (b) Pc; (c) MSSOR.

Fig. 10. Eigenvalue distribution of preconditioned matrices using the largest finite element mesh with 2020 elements for the plane strain pile group problem: (a) GJ preconditioned matrix; (b) Pc preconditioned matrix; (c) MSSOR preconditioned matrix.

Table 2 Convergence behaviors of MSSOR preconditioned SQMR method based on different criteria with the tolerance 1E6. Nels Case b1, iters/R 661 1116 2020 Case b2, iters/R 661 1116 2020 Case b3, iters/R 661 1116 2020

m + n (N)

Rr

Ri

Rus & Rps

Rui & Rpi

3980 + 684 6716 + 1148 12124 + 2069

874/1.5E7 1172/9.8E8 2124/1.0E7

174/9.3E1 235/9.4E1 1173/5.3E1

873/2.2E7 1147/2.7E7 2122/1.2E7

624/9.5E2 235/9.4E1 1173/5.3E1

3980 + 684 6716 + 1148 12124 + 2069

1065/7.0E9 1448/8.0E9 2923/1.9E8

307/1.3E1 477/1.5E1 737/3.9E1

877/2.7E6 1282/6.0E7 2437/1.6E6

307/1.3E1 477/1.5E1 737/3.9E1

3980 + 684 6716 + 1148 12124 + 2069

858/1.7E7 1178/1.6E7 2246/1.4E7

768/1.4E5 848/1.5E2 1455/8.5E2

861/1.2E7 1178/1.6E7 2246/1.4E7

768/1.4E5 848/1.5E2 1455/8.5E2

1281

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

several orders higher than the relative errors due to the relative residual norm criteria, indicating that Ri is less reliable. Further, it can be seen that apart from Ri falling below R, the fairly sizeable local oscillations exacerbate this problem. By comparing the sub-plots in Figs. 7–9, it can be observed that Pc preconditioned SQMR has comparatively smoother convergence behaviors compared to GJ and MSSOR preconditioned SQMR. This difference may be adequately interpreted by the spectral properties of these preconditioned matrices. As shown in Fig. 10 for eigenvalue distributions, both GJ and MSSOR preconditioned matrices have imaginary parts in their eigenvalues, but the Pc preconditioned matrix only has real eigenvalues. We also notice that in Fig. 6 the eigenvalues for the tensile plate example are real and the convergence curves in Figs. 4 and 5 are relatively smooth.

Hence, it is postulated here that the presence of imaginary parts in the eigenvalues may create or amplify the oscillatory convergence, which is undesirable from a stopping criterion point of view. The imaginary parts for GJ are larger than those for MSSOR. This may be correlated to the slightly more severe oscillations in Figs. 7a, 8a, and 9a. Fig. 11 examines the effect of initial guess x0 on the Pc preconditioned SQMR method. It can be seen that for the zero initial guess, i.e. x0 = zeros(N, 1), the Rr convergence curve lies above the Ri convergence curve. For two non-zero initial guesses, the Rr con-

104 10

10

4

10

1

R Ri

10

-2

Rr

10

-5

||rk||2

10-8

10

-2

Rr

10

-5

||rk||2

10

10

-14

-8

-11

10

-17

10 10

-14

10

-17

(a) 0

500

104 10

1000 1500 2000 2500 R Ri

1

10

-2

Rr

10

-5

||rk||2

10

-17

10

-20

(b) 0 10

4

10

1

10

4

10

1

3000 R Ri

10-2

Rr

10-5

||rk||2

10

-8

10-17 10-20 (b) 0 500

1000 1500 2000 2500 R Ri

-2

Rr ||rk||2

10-5 10-8

2000

4000

104

6000

1

R Ri

10-2

Rr

10

10

-5

10

-8

||rk||2

10-11

10-11

10-14

10-14

-17

10

-17

10

10

-20

10-20 (c) 0

(c) 0

2000

10-14

10-11 -14

10

1000

10-11

10-8

10

10-11

10-20 (a) 0

10-20

R Ri

1

500

1000 1500 2000 2500

Fig. 11. Convergence behaviors of SQMR preconditioned by Pc preconditioner based on global criteria for the plane strain pile group problem (b3 RHS case): (a) x0 = zeros(N, 1); (b) x0 = ones(N, 1); (c) x0 = rand(N, 1).

2000

4000

6000

Fig. 12. Convergence behaviors of SQMR preconditioned by MSSOR preconditioner based on global criteria for the plane strain pile group problem (b3 RHS case): (a) x0 = zeros(N, 1); (b) x0 = ones(N, 1); (c) x0 = rand(N, 1).

1282

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

vergence curve lies below Ri the convergence curve. Fig. 12 examines the effect of initial guess x0 on the MSSOR preconditioned SQMR method. When the SQMR solver is terminated at Ri 6 1E6, the required iteration is 1455 if x0 = zeros(N, 1), and the corresponding R and ||rk||2 are 8.48E2 and 3.56E+1, respectively. The required iteration is 1613 if x0 = ones(N, 1), and the corresponding R and ||rk||2 are 2.52E+1 and 1.45E+4, respectively. The required iteration is 2210 if x0 = rand(N, 1), and the corresponding R and ||rk||2 are 3.47E4 and 1.24E1, respectively. When the SQMR solver is terminated at Rr 6 1E6, the required iteration is 2246 if x0 = zeros(N, 1), and the corresponding R and ||rk||2 are 1.43E7 and 3.70E7, respectively. The required iteration is 2127 if x0 = ones(N, 1), and the corresponding R and ||rk||2 are 8.89E1 and 3.58E1, respectively. The required iteration is

103 0

R Rsu

10

-3

Rsp

10

-6

10

10

R

u r

Rrp

-9

10

-12

10

-15

10

-18

2043 if x0 = rand(N, 1), and the corresponding R and ||rk||2 are 1.70E+0 and 1.09E+0, respectively. From Figs. 11 and 12, we know that for the zero initial guess (usually the default choice), Rr is stricter than Ri. Further, it can be observed that Ri should be carefully evaluated if the eigenvalues of preconditioned matrix are complex, because the dramatic local oscillations may lead to premature termination of an iterative solver. 5.2.2. Separated and decoupled convergence criteria Unlike the initial guess, RHS does not have an obvious influence on the convergence behaviors of an iterative solver as shown in Figs. 7–9. Hence, this section only presents results for the b3 case. Fig. 13a and b show that the separated and the decoupled criteria basically follow the convergence behaviors of their global counterparts (Rr and Ri in Fig. 9c). Because the convergence characteristics between the global and the separated criteria are comparable, it is computationally advantageous to check one of the separated

10

-1

10

-4

10

-7

R Ri Rr

10

-10

10

-13

||rk||2

10-16

(a) 0

1000 2000 3000 4000 5000

0

R Riu

10

-3

Rip

10

-6

10

-9

10

10-18

103

10

-4

10

-7

1000 2000 3000 4000 5000

10

-16

10

-19

10

-22

Rr ||rk||2

(b) 0 10

-1

10

-3

Rp

10

-4

10

-7

10

-6

10

-9

10-13

-15

10

-16

10 10

10

-19

-18

Fig. 13. Convergence behaviors of MSSOR preconditioned SQMR based on the separated or decoupled criteria for the plane strain pile group problem (b3 RHS case).

1000 R Ri

||rk||2

-10

-12

1000 2000 3000 4000 5000

500

Rr

10

(c) 0

R Ri

-10

0

10

1000

10-1

R Ru

10

500

10-13

-15

(b) 0

(a) 0

10

10-12 10

-19

10-22

3

10

10

10-22

(c) 0

50

100

150

200

250

Fig. 14. Convergence behaviors of preconditioned iterative solvers based on different criteria for the 3-D footing problem: (a) JCG; (b) JSQMR; (c) ICCG(0).

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

criteria only (i.e. Rui or Rpi ; Rur or Rpr ). For the decoupled criteria u (Rus and Rps Þ; theev aluationsofRs and Rps need additional computations, and thus checking the decoupled criteria every nth iterations should be more practical. There are two additional points of interest. Mathematically, Rur is equal to Rus . However, their convergence profiles in floating-point (finite precision) arithmetic are distinct. The decoupled criteria produce smooth convergence profiles, which is also advantageous as noted previously. We remarked in Section 4 that the separated/decoupled residual vectors are shorter than the global vector. In particular, the excess pore water pressure vector may only constitute 10% of the global vector. Fig. 13c compares the separated relative error norm, Ru and Rp, with the global ‘‘ideal” relative error norm R. Despite Rp being computed from a significantly shorter vector, it tracks R even more closely than Ru. Overall, it would appear that the 2-norm is

10

-1

10

-4

10

-7

10

R Ri Rr ||rk||2

-10

10-13 10-16 10-19 10-22

(a) 0 10

-1

10

-4

10

-7

1000

10 10

-13

Rr

10-19 -22

(b) 0

500

10-4 10

1000 R Ri

10-1

Rr

-7

||rk||2

10-10 10

-13

10-16 10-19 10-22

(c) 0

500

5.3. A 3-D footing foundation example A 3-D footing foundation problem as shown in Fig. 3 is analyzed with 1440 elements under drained and consolidation conditions. A symmetric positive definite system (N = 17,980) and a symmetric indefinite system (N = 19,670) are generated, respectively. The heterogeneous soil body is simulated with E’ = 50 MPa, v0 = 0.3 for upper soil layer and E0 = 1 MPa, v0 = 0.3 for lower soil layer. For consolidation analysis, k = 1E5 m/s is simulated for upper soil layer and k = 1E9 m/s is simulated for lower soil layer. Fig. 14. shows the convergence behaviors of different preconditioned iterative solvers for a zero initial guess vector [i.e. x0 = zeros(N, 1)]. It can be observed that the convergence curve due to ||rk||2 is almost same as that due to Rr, indicating that in this case ||r0||2 = ||b||2 1. In addition, for this symmetric positive definite system, the convergence behaviors of an iterative solver based on Rr and Ri are similar and the convergence curves are relatively smooth because the eigenvalues of these preconditioned matrices are all real. Fig. 15 provides the convergence behaviors of SQMR preconditioned by GJ, Pc and MSSOR, respectively. In contrast to the pile-group consolidation analysis (Fig. 7), the Rr convergence curve lies below the R curve. The observations that R may lie above or below Rr are theoretically supported by Eq. (11). For the zero initial guess which is usually a default selection, the convergence curve due to Rr criterion lies above that due to Ri criterion as illustrated by Figs. 7–9 and 15. Hence, for such an initial guess, Rr criterion should be stricter than Ri.

6. Conclusion

R Ri

10-16

10

robust against gross differences in vector length and it is feasible to use either one of the separated/decoupled criteria as a stopping criterion.

2000

||rk||2

-10

1283

1000

Fig. 15. Convergence behaviors of SQMR preconditioned by different preconditioners based on global criteria for the 3-D footing problem (b1 RHS case): (a) GJ; (b) Pc; (c) MSSOR.

In this study, several popular convergence criteria are evaluated for two types of linear systems of equations: symmetric positive definite linear system and symmetric indefinite linear system. The former is illustrated using a tensile plate problem and drained analyses of a 2-D pile group and a 3-D footing. The latter is illustrated using coupled consolidation analyses of a 2-D pile group and a 3-D footing. We can summarize the following observations and remarks: (1) When the penalty method is applied for the constrained DOFs, the standard relative ‘‘improvement” norm criterion Ri and the relative residual norm criterion Rr will result in premature termination. The recommended solution is to remove the components corresponding to the restrained DOFs. (2) The initial guess may influence the convergence behaviors of an iterative solver. The difference between R and Rr (or Ri) changes with the initial guess vector. For the commonly adopted zero initial guess vector, a symmetric iterative solver such as PCG or SQMR exhibits similar convergence profiles based on Ri and Rr for the symmetric positive definite linear system. For the symmetric indefinite system, the convergence curves based on Rr and Ri are apparently more affected by the initial guess. (3) For the symmetric indefinite linear system with zero initial guess, Ri tends to produce larger error when referenced to R for a prescribed tolerance. The reverse is true for Rr. A tolerance of 1E6 appears to be reasonable for Rr. However, a tolerance significantly smaller than 1E6 should be used in practical finite element computations if Ri is selected as the stopping criterion.

1284

X. Chen, K.K. Phoon / Computers and Geotechnics 36 (2009) 1272–1284

(4) For the symmetric indefinite linear system, the convergence profiles may exhibit dramatic local oscillations. If one of these oscillations dips below the prescribed tolerance momentarily during the initial stage of the iterative process, a fairly extreme premature termination will take place for iterative solvers based on Ri criterion. Numerical evidence appears to show that these dramatic oscillations are present when the eigenvalues of the preconditioned matrix carry imaginary components. (5) For the consolidation problem with two types of variables, the separated or the decoupled criteria can track the global criteria very well, indicating that checking one of the separated or the decoupled criteria is a practical and cheaper option. It is interesting that the decoupled relative residual norm criteria track the relative error norm even more closely than the global relative residual norm. They have the added advantage of exhibiting smooth profiles. However, the evaluations of the decoupled criteria need additional computations, and thus checking these criteria every nth iterations should be more practical.

References [1] Axelsson O, Kaporin I. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. Numer Linear Algebra Appl 2001;8(4):265–86. [2] Arioli M. Stopping criteria for iterations in finite element methods. Numer Math 2005;99(3):381–410. [3] Arioli M, Loghin D. Stopping criteria for mixed finite element problems. RAL Technical Reports 2006, Technical Report, RAL-TR-2006-010. [4] Smith IM, Griffiths DV. Programming the finite element method. 4th ed. John Wiley & Sons; 2004. [5] Picasso M. A stopping criterion for the conjugate gradient algorithm in the framework of anisotropic adaptive finite elements. Commun Numer Methods Eng 2008 [Early View]. [6] Barrett R, Berry M, Chan T, Demmel J, Donato J, Dongarra J, et al. Templates for the solution of linear systems: building blocks for iterative methods. Philadelphia, PA: SIAM Press; 1994. [7] Chen X. Preconditioners for iterative solutions of large-scale linear systems arising from Biot’s consolidation equations. Ph.D. Thesis, Department of Civil Engineering, National University of Singapore; 2005. [8] Driscoll TA, Toh KC, Trefethen LN. From potential theory to matrix iterations in six steps. SIAM Rev 1998;40:547–78. [9] Greenbaum A. Iterative methods for solving linear systems. Philadelphia: SIAM; 1997. [10] Trefethen LN, Embree M. Spectra and pseudospectra: the behavior of non-normal matrices and operators. New Jersey: Princeton University Press; 2005. [11] Saad Y. Iterative methods for sparse linear systems. Boston: PWS Publishing Company; 1996. [12] van der Vorst HA. Iterative Krylov methods for large linear systems. New York: Cambridge University Press; 2003. [13] Toh KC. GMRES vs. ideal GMRES. SIAM J Matrix Anal Appl 1997;18(1):30–6. [14] Simoncini V, Szyld DB. Recent computational developments in Krylov subspace methods for linear systems. Numer Linear Algebra Appl 2007; 14(1):1–59. [15] Nachtigal NM, Reddy SC, Trefethen LN. How fast are nonsymmetric matrix iterations? SIAM J Matrix Anal Appl 1992;13(3):778–95.

[16] Freund RW, Nachtigal NM. A new Krylov-subspace method for symmetric indefinite linear systems. In: Ames WF, editor. Proceedings of the 14th IMACS world congress on computation and applied mathematics. IMACS; 1994. p. 1253–6. [17] Kincaid DR, Respess JR, Young DM, Grimes RG. Algorithm 586: ITPACK 2C: a FORTRAN package for solving large sparse linear systems by adaptive accelerated iterative methods. ACM Trans Math Softw 1982;8(3):302–22. [18] Freund RW. Preconditioning of symmetric, but highly indefinite linear systems. In: Proceedings of the 15th IMACS world congress on scientific computation, modelling and applied mathematics, vol. 2; 1997. p. 551–6. [19] Benzi M. Preconditioning techniques for large linear systems: a survey. J Comput Phys 2002;182:418–77. [20] Chen K. Matrix preconditioning techniques and applications. Cambridge monographs in applied and computational mathematics, vol. 19. Cambridge: Cambridge University Press; 2005. [21] Mroueh H, Shahrour I. Use of sparse iterative methods for the resolution of three-dimensional soil/structure interaction problems. Int J Numer Anal Methods Geomech 1999;23(15):1961–75. [22] Britto AM, Gunn MJ. Critical state soil mechanics via finite elements. Chichester, West Sussex: Ellis Horwood Ltd.; 1987. [23] Biot MA. General theory of three-dimensional consolidation. J Appl Phys 1941;12:155–64. [24] Onate E, Rojek J, Taylor RL, Zienkiewicz OC. Finite calculus formulation for incompressible solids using linear triangles and tetrahedra. Int J Numer Methods Eng 2004;59(11):1473–500. [25] Smith IM. A general purpose system for finite element analyses in parallel. Eng Comput 2000;17(1):75–91. [26] Chen X, Phoon KK, Toh KC. Partitioned versus global Krylov subspace iterative methods for FE solution of 3-D Biot’s problem. Comput Methods Appl Mech Eng 2007;196(25–28):2737–50. [27] Bergamaschia L, Ferronato M, Gambolati G. Novel preconditioners for the iterative solution to FE-discretized coupled consolidation equations. Comput Methods Appl Mech Eng 2007;196(25–28):2647–56. [28] Chen X, Phoon KK, Toh KC. Performance of zero-level fill-in preconditioning techniques for iterative solutions in geotechnical applications, submitted for publication. [29] GeoSoft, GeoFEA Users’ Manual 8.0, Singapore; 2006. [30] Toh KC, Phoon KK. Comparison between iterative solution of symmetric and non-symmetric forms of Biot’s FEM equations using the generalized Jacobi preconditioner. Int J Numer Anal Methods Geomech 2008;32(9):1131–46. [31] Ferronato M, Pini G, Gambolati G. The role of preconditioning in the solution to FE coupled consolidation equations by Krylov subspace methods. Int J Numer Anal Methods Geomech 2008;33(3):405–23. [32] Chen X, Toh KC, Phoon KK. A modified SSOR preconditioner for sparse symmetric indefinite linear systems of equations. Int J Numer Methods Eng 2006;65(6):785–807. [33] Phoon KK, Toh KC, Chan SH, Lee FH. An efficient diagonal preconditioner for finite element solution of Biot’s consolidation equations. Int J Numer Methods Eng 2002;55(4):377–400. [34] Elman HC, Silvester DJ, Wathen AJ. Iterative methods for problems in computational fluid dynamics. Technical Report 96/19, Oxford University Computing Laboratory; 1996. [35] van der Vorst HA. Efficient and reliable iterative methods for linear systems. J Comput Appl Math 2002;149(1):251–65. [36] Meijerink JA, van der Vorst HA. An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix. Math Comput 1977;31(137):148–62. [37] Zhang Y. Solving large-scale linear programs by interior-point methods under the Matlab environment. Optim Methods Softw 1998;10(1):1–31. [38] Krabbenhoft K, Lyamin AV, Sloan SW, Wriggers P. An interior-point algorithm for elastoplasticity. Int J Numer Methods Eng 2007;69(3):592–626. [39] Axelsson O. Iteration number for the conjugate gradient method. Math Comput Simulat 2003;61(3–6):421–35. [40] Phoon KK, Toh KC, Chen X. Block constrained versus generalized Jacobi preconditioners for iterative solution of large-scale Biot’s FEM equations. Comput Struct 2004;82(28):2401–11.