A domain decomposition method of stochastic PDEs: An iterative solution techniques using a two-level scalable preconditioner

Journal of Computational Physics 257 (2014) 298–317

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

A domain decomposition method of stochastic PDEs: An iterative solution techniques using a two-level scalable preconditioner ✩ Waad Subber, Abhijit Sarkar ∗ Department of Civil and Environmental Engineering, Carleton University, Ottawa, Ontario K1S5B6, Canada

a r t i c l e

i n f o

Article history: Received 28 July 2012 Received in revised form 16 August 2013 Accepted 29 August 2013 Available online 7 September 2013 Keywords: Domain decomposition method Schur complement system Neumann–Neumann preconditioner Balancing domain decomposition by constraints Stochastic PDEs Polynomial chaos expansion Stochastic finite element method

a b s t r a c t Recent advances in high performance computing systems and sensing technologies motivate computational simulations with extremely high resolution models with capabilities to quantify uncertainties for credible numerical predictions. A two-level domain decomposition method is reported in this investigation to devise a linear solver for the large-scale system in the Galerkin spectral stochastic finite element method (SSFEM). In particular, a two-level scalable preconditioner is introduced in order to iteratively solve the large-scale linear system in the intrusive SSFEM using an iterative substructuring based domain decomposition solver. The implementation of the algorithm involves solving a local problem on each subdomain that constructs the local part of the preconditioner and a coarse problem that propagates information globally among the subdomains. The numerical and parallel scalabilities of the two-level preconditioner are contrasted with the previously developed one-level preconditioner for two-dimensional flow through porous media and elasticity problems with spatially varying non-Gaussian material properties. A distributed implementation of the parallel algorithm is carried out using MPI and PETSc parallel libraries. The scalabilities of the algorithm are investigated in a Linux cluster. © 2013 Elsevier Inc. All rights reserved.

1. Introduction In the numerical modeling of many practical engineering and natural systems, the random heterogeneities and imperfections of the propagating media may significantly influence the predictive capabilities of computer simulations. Recent advances of high performance computing and sensing technologies motivate computational simulations with extremely high resolution which should integrate efficient uncertainty quantification methods for realistic and credible numerical predictions. For such extreme scale simulations, uncertainty quantification using the standard Monte Carlo simulation may be time-consuming or impractical. To reduce the computational cost of the traditional Monte Carlo simulation, a multi-level Monte Carlo method has recently been proposed [2,3]. Alternatives to Monte Carlo sampling techniques are the spectral stochastic finite element method (SSFEM) whose implementations can exploit non-intrusive collocation method (e.g. [4–11]) or intrusive Galerkin approach (e.g. [12–21]). This investigation focuses on an efficient distributed implementation of Galerkin scheme in SSFEM for large-scale systems, circumventing the need of any random or deterministic sampling. For a comprehensive review of SSFEM, we refer to e.g. [12,14,22–26].

✩

*

The preliminary version of the paper is published in the proceeding of HPCS 2011 conference [1]. Corresponding author. E-mail addresses: [email protected] (W. Subber), [email protected] (A. Sarkar).

0021-9991/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.jcp.2013.08.058

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

299

In the so-called intrusive SSFEM (e.g. [12]), a Galerkin projection scheme transforms the original stochastic linear system into a coupled set of block linear systems. Depending on the scales of random fluctuations present on the system parameters and degree of nonlinearities between the input and output processes, the order of the resulting deterministic linear system can increase significantly in the SSFEM. The computational efficiency of the SSFEM solver is primarily dictated by the solution strategies adopted to tackle the coupled deterministic linear systems. This fact inspires the development of a scalable solver for the SSFEM which can efficiently exploit the modern high performance computing systems. Preconditioned iterative solvers [27–32] and multigrid methods [33–35,26] have been exploited to develop robust and reliable iterative solvers for the Galerkin projection-based intrusive SSFEM. The block-Jacobi preconditioner based on the mean stiffness matrix is generally employed to accelerate the convergence rate of the iterative methods. However, such iterative solvers equipped with the block-Jacobi preconditioner show performance degradation as the level of uncertainty and order of the deterministic linear system grow in the SSFEM [27,28,30,32,36]. This fact prompts further research to enhance the robustness and performance of iterative techniques in order to solve the large-scale system in the intrusive SSFEM. To this end, a non-overlapping domain decomposition algorithm is reported in [37] using the intrusive SSFEM. The algorithm relies on a Schur complement-based domain decomposition in the physical space and a polynomial expansion of stochastic processes representing the system parameters and solution field. Extending the basic mathematical formulation in [37], a one-level iterative substructuring technique is presented in [38,39] for efficient solution of the large-scale linear system in the intrusive SSFEM. Numerical experiments demonstrate that the convergence rate of the iterative algorithm is mainly insensitive to the magnitude of the coefficient of variation (i.e. strength of randomness) of the system parameters and order of stochastic dimension. However, the one-level iterative algorithm in [38,39] shows that the iteration count grows linearly as we increase the number of subdomain partitions. A two-level domain decomposition preconditioner is described in this paper to improve the efficiency of the iterative substructuring technique in the intrusive SSFEM [38,39]. A dual–primal domain decomposition method for SSFEM is reported elsewhere by the authors [40–43]. In particular, the one-level Neumann–Neumann (NN) preconditioner in [38,39] is complemented by a coarse grid in this paper. A collection of corner nodes at the subdomain interfaces constitute the coarse grid. Consequently, a coarse problem is solved at each iteration to spread information across all the subdomains. This information exchange achieved by the coarse grid makes the preconditioner scalable. The two-level preconditioner may be construed to be a probabilistic extension of the Balancing Domain Decomposition by Constraints (BDDC) [44,45] devised for deterministic PDEs. The parallel performances of the previously developed one-level preconditioner [38,39] and the new two-level preconditioner are contrasted using illustrations from linear elasticity and flow through porous media with spatially varying non-Gaussian material properties. We use PETSc [46] and MPI [47] parallel libraries for the distributed implementation of the algorithm. 2. Schur complement system of the stochastic PDEs We briefly describe Schur complement-based domain decomposition solver for stochastic PDEs [1,37–39,48,40,41]. For an elementary exposition of the methodology, we consider an elliptic stochastic PDE on a domain Ω with a prescribed homogeneous Dirichlet boundary condition on ∂Ω as

  ∇ · κ (x, θ)∇ u (x, θ) = f (x), u (x, θ) = 0,

x ∈ Ω,

x ∈ ∂Ω,

(1) (2)

where θ denotes an element in the sample space, u (x, θ) is the solution process, κ (x, θ) represents the diffusion coefficient modeled as a strictly positive random field and f (x) is the deterministic forcing term. The finite element discretization of the above elliptic stochastic PDE yields the following stochastic linear system

A(θ)u(θ) = f,

(3)

where A(θ) is the stochastic stiffness matrix, u(θ) is the random response vector and f is the external forcing vector. For high resolution models, domain decomposition techniques can be exploited to tackle Eq. (3) on parallel computers [1,37–39, 48,40,41]. ns The computational domain Ω is partitioned into ns non-overlapping subdomains Ω = s= 1 Ωs with interfaces defined ns as Γ = s=1 Γs where Γs = ∂Ωs \∂Ω [49–52]. The nodes of each subdomain (Ωs , s = 1, . . . , ns ) are partitioned into two subsets: interior nodes that belong to the given subdomain and interface nodes shared by two or more adjacent subdomains as shown in Fig. 1. According to this partition, the equilibrium system of a typical subdomain Ωs is expressed as



AIIs (θ) s AΓ I (θ)

AsIΓ (θ)

s AΓΓ (θ)



usI (θ) s uΓ (θ)





=

fsI s fΓ



(4)

s where usI (θ) denotes the unknowns associated with the interior nodes of the subdomain Ωs and uΓ (θ) corresponds to the interface unknowns associated with the interface boundary Γs .

300

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 1. Mesh decomposition with interior nodes ( ) and interface nodes ( ).

The non-Gaussian random material properties of the system are represented by the polynomial chaos expansion (PCE) [12,25,22]. The subdomain equilibrium equation can then be expressed as L



AIIs ,i

Ψi (θ)

AsIΓ,i

s AΓ I ,i

i =0



s AΓΓ ,i

usI (θ)



 =

s uΓ (θ)

fsI s fΓ

 (5)

,

where Ψi (θ) is a set of multi-dimensional Hermite polynomials with orthogonality property defined as





Ψi (θ)Ψ j (θ) = Ψi2 (θ) δi j ,

(6)

where δi j defines the Kronecker delta and · represents the expectation operator. We expand the solution vector of the subdomain using the PCE as



usI (θ)

 =

s uΓ (θ)

N

 Ψ j (θ)

j =0

usI, j

s uΓ, j

 (7)

.

Performing a global assembly of the interface unknowns and then performing Galerkin projection on the PC basis, we obtain the following linear system (see [1,37–39,48,40,41] for comprehensive details)

⎡

A1II .. .

⎢ ⎢ ⎢ ⎢ ⎣

... .. .

0

.. . n ... AIIs . . . RnTs AnΓsI

0

R1T A1Γ I

⎤⎧ 1 ⎪ ⎪ UI ⎥⎪ ⎪ ⎥ ⎨ .. . ⎥ ⎥ ⎪ ns n ⎦⎪ A I Γs Rns U ⎪ ⎪ ns ⎩ I T s UΓ R A R s s ΓΓ s =1 A1IΓ R1 .. .

⎫ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎭

=

⎧ ⎪ ⎪ ⎪ ⎪ ⎨

F I1 .. .

⎫ ⎪ ⎪ ⎪ ⎪ ⎬

n ⎪ ⎪ ⎪ ⎪ F s ⎪ ⎪ ⎪ ⎭ ⎩ ns I T s ⎪ R F s Γ s =1

,

(8)

where L  s  Aα β jk = Ψi Ψ j Ψk Aαs β,i ,

U Im

=

i =0

T m um , I ,0 , . . . , u I , N



The subscripts

Fαs ,k = Ψk fαs ,

UΓ = (uΓ,0 , . . . , uΓ, N )T .

α and β denote the index I and Γ and the restriction operator Rs is defined as

Rs = blockdiag(Rs , . . . , Rs ), !" #

(9)

N +1

where N + 1 is the number of terms in the PCE and Rs is a Boolean rectangular matrix representing a scatter operator of s (θ) as the global interface vector uΓ (θ) into its local components uΓ s uΓ (θ) = Rs uΓ (θ).

(10)

Eliminating the interior variables in Eq. (8), the global extended Schur complement system becomes

SUΓ = GΓ ,

(11)

where

S=

ns

 s  − 1 s  RsT AΓΓ − AΓs I AIIs A IΓ Rs ,

(12)

s =1

GΓ =

ns

  − 1 s  RsT FΓs − AΓs I AIIs FI .

(13)

s =1

The extended Schur complement system in Eq. (11) is solved iteratively and in parallel to obtain the interface variables UΓ . These variables relate to the interface nodes as schematized in Fig. 2.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

301

Fig. 2. The interface nodes ( ) for each PCE coefficient of the solution process (e.g., 3 coefficients).

3. Iterative solution of the stochastic interface system In the iterative substructuring technique, the interface system given by Eq. (11) is generally tackled by the preconditioned conjugate gradient method (PCGM). In the PCGM, a preconditioning matrix M−1 is employed to transform the original system into an equivalent system whereby the transformed system possesses a smaller condition number than that of the original system. Additionally the preconditioned system has eigenvalues clustered near one [53–55]. The preconditioned extended Schur complement system is

M−1 SUΓ = M−1 GΓ .

(14)

In this paper, we devise a two-level preconditioner to enhance the performance of iterative substructuring techniques to tackle Eq. (11). The preconditioner consists of a local part to smooth some components of the error and a global part to propagate the error globally over the subdomains. The general form of the two-level preconditioner is

M−1 =

ns

−1 1 T RsT M− s Rs + R0 M0 R0 .

s =1

Within the parallel PCGM algorithm, the two-level preconditioner can be implemented using Algorithm 1 [38,40]. Algorithm 1: The parallel PCGM procedure.

5

: (UΓ0 ) s 0 T compute : r0Γ = GΓ − ns= 1 R s S s R s UΓ  ns 0 T − 1 precondition : Z = [ s=1 Rs Ms Rs + R0T M0−1 R0 ]r0Γ compute : P0 = Z0 compute : ρ 0 = (r0Γ , Z 0 ) for j = 0, 1, 2, . . . , until convergence:

6

do

input 1 2 3 4

7

compute

: Qj = j tmp j

ns

T j s=1 Rs Ss Rs P j j

= (Q , P )

8

compute

:ρ

9

compute

: α = ρ j /ρtmp

update

: UΓ

10 11 12 13 14 15

j +1

j +1

j

j

= UΓ + α j P j j

= rΓ − α j Q j  j +1 −1 T 1 precondition : Z j +1 = [ ns=s 1 RsT M− s Rs + R0 M0 R0 ]rΓ update

: rΓ

compute compute update

: ρ j +1 = (rΓ , Z j +1 ) : β j = ρ j +1 /ρ j : P j +1 = Z j +1 + β j P j

output

j +1

: (UΓ )

(15)

302

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Neither the extended Schur complement matrix in step 7, nor the two-level preconditioner in step 12 need to be formed explicitly. The effect of the extended Schur complement matrix on a given vector is obtained by computing three subdomain level matrix vector products and solving local stochastic Dirichlet problems as outlined in Algorithm 2 [38,40]. Whereas the preconditioner effect on the residual is calculated by solving local stochastic Neumann problems and a global coarse problem. Algorithm 2: Dirichlet-solver procedure.

1 2

input : (P ) for s = 1, 2, . . . , ns in parallel: do

scatter : P s = Rs P compute : v 1 = AsIΓ P s solve : AIIs v 2 = v 1 compute : v 3 = AΓs I v 2 s compute : v 4 = AΓΓ Ps s compute : Q = v 4 − v 3 s T s gather : Q = ns= 1 Rs Q

3 4 5 6 7 8 9

output

: (Q)

4. Preconditioners for stochastic PDEs For stochastic PDEs, the convergence rate of the PCGM equipped with efficient preconditioners should not be hampered by the level of uncertainty and order of the stochastic dimension. Furthermore, the PCG solvers must also demonstrate a scalable performance in relation to the physical problem size and number of subdomains. Such performance can be achieved by a two-level preconditioner suitably customized for the stochastic system as delineated next. 4.1. One-level preconditioner Previously, the authors introduced a one-level Neumann–Neumann preconditioner for stochastic PDEs [38,39]. The extended Neumann–Neumann preconditioner is defined as 1 M− NN =

ns

RsT DsT Ss−1 Ds Rs ,

(16)

s =1

where

Ds = blockdiag(Ds , . . . , Ds ), !" #

(17)

N +1

and Ds is a diagonal scaling matrix defined by ns

RsT Ds Rs = I.

(18)

s =1

The diagonal elements of Ds are equal to the reciprocal of the number of subdomains having the same interface nodes [49–52]. The implementation of the extended Neumann–Neumann preconditioner entails solving a subdomain level stochastic Neumann problem. This step calculates the effect of the inverse of the local Schur complement matrix Ss−1 on the local s residual vector rΓ as



AIIs AΓs I



AsIΓ s AΓΓ j

Xs UΓs





= j

0 s rΓ



,

(19)

s where rΓ = Ds Rs rΓ and rΓ denotes the global residual at the jth iteration of the PCGM. The PCGM equipped with a one-level Neumann–Neumann preconditioner exhibits a convergence rate that is almost insensitive to the level of uncertainty in the system parameters [38,39]. As no mechanism for global information exchange exists in the one-level Neumann–Neumann preconditioner, its performance deteriorates for relatively large number of subdomains. For scalability, the global communication is required to propagate information across the subdomains [49–52]. Another difficulty in the implementation of the one-level Neumann–Neumann preconditioner is that the local stochastic

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

303

Fig. 3. Partitioned interface nodes of each PCE coefficient of the solution process (e.g., 3 coefficients) into: remaining ( ) and corner ( ) nodes.

Neumann problem defined in Eq. (19) is singular for floating subdomains. In practice, the singularity of the local stiffness matrix is avoided by altering its diagonal entries with a small positive constant and the local Neumann problem is solved approximately [51,52]. However, such approximation may adversely affect the overall convergence rate of the algorithm [51]. 4.2. Two-level preconditioner In the one-level Neumann–Neumann preconditioner, the information is only exchanged among the neighboring subdomains. To provide a mechanism for a global transfer of information, a coarse problem is required [49,50,52]. To this end, we describe a mathematical framework for a two-level preconditioner for the iterative substructuring to devise an iterative solver for the large-scale linear system in the intrusive SSFEM. This method bears similarity with a version of the Balancing Domain Decomposition by Constraints (BDDC) for deterministic PDEs [44,45]. The two-level preconditioner is formulated by adding a coarse grid to the one-level Neumann–Neumann preconditioner as described next. ns The interface nodes of each subdomain {Ωs }s= 1 are partitioned into a set of remaining nodes (boundary nodes shared only between two adjacent subdomains) and corner nodes (boundary nodes shared among more than two subdomains plus the nodes at the ends of interface edges) [44,45,56,57] as shown in Fig. 3. Let Rrs and Rcs be the Boolean operators that split the local interface unknowns UΓs into remaining Urs and corner Ucs variables respectively as



Urs Ucs





=

Rrs Rcs



UΓs

(20)

where

  Rrs = blockdiag Rrs , . . . , Rrs , !" # Rcs

= blockdiag



N +1

Rcs , . . . , Rcs

!"

N +1

#



.

For the kth polynomial chaos coefficients, the matrices Rrs and Rcs are s urs,k = Rrs uΓ, k, s ucs ,k = Rcs uΓ, k.

According to the partition defined in Eq. (20), the subdomain level Neumann problem given by Eq. (19), can be rewritten as

⎡

Aiis ⎣ As ri s Aci

s Air s Arr s Acr

where

Frs = Rrs rΓs , Fcs = Rcs rΓs .

⎧ ⎫ s ⎤⎧ s ⎫ Aic ⎨X ⎬ ⎨ 0 ⎬ s ⎦ Arc Us = Fs , ⎩ rs ⎭ ⎩ rs ⎭ s Acc Uc Fc

(21)

304

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Performing a partial assembly at the corner nodes for continuity, leads to the following linear system

⎡ ⎢ ⎣

ns

Aiis s Ari

sT s s=1 Bc Aci

s Air s Arr

ns

sT s s=1 Bc Acr

⎧ ⎫ ⎤⎧ s ⎫ 0 ⎪ ⎨ ⎬ ⎨X ⎬ ⎪ ⎥ s s Fr , ⎦ Ur = ⎭ ⎪ ⎩ ns ⎩ ns B s T F s ⎪ ⎭ sT s s U B A B c cc c c s =1 c s =1 c s s Aic Bc s Arc Bcs

(22)

n

s where Bcs is a restriction operator that maps the global corner unknown Uc into its local components {Ucs }s= 1 as

Ucs = Bcs Uc ,

(23)

where

  Bcs = blockdiag Bcs , . . . , Bcs , !" # N +1

and Bcs is defined by

ucs ,k = Bcs uc ,k . Eliminating the local variables X s and Urs in Eq. (22), we obtain the following symmetric positive definite coarse problem

F cc Uc = dc ,

(24)

where

F cc =

ns

Bcs

T s Scc

 s −1 s  s s − Scr Srr Src Bc ,

s =1

dc =

ns

Bcs

T

 s −1 s  s Fcs − Scr Srr Fr ,

s =1 s s s −1 s and Sαs β = Aα β − Aα i [Aii ] Ai β where α and β stand for the index r and c. The corner nodes artificially enforce Dirichlet boundary conditions to eliminate the singularity arising from the subdos main Schur complement matrix Srr of the floating subdomains. The remaining interface variables Urs can be obtained next for each subdomain as

s Srrs Urs = Frs − Src Ucs .

(25)

Once the local interface variables Urs and Ucs are available, the global interface unknowns UΓ can be aggregated as

UΓ =

ns

 T  T RsT Ds Rrs Urs + Rcs Ucs .

(26)

s =1

Substituting Eq. (24) and Eq. (25) into Eq. (26) and performing some algebraic manipulations, the two-level preconditioner is written as 1 M− NNC =

ns

RsT Ds Rrs

T  s  −1 r Srr Rs Ds Rs

+ R 0T [ F cc ]−1 R 0 ,

(27)

s =1

where the coarse scale restriction matrix R 0 takes the following form

R0 =

ns

Bcs

T

 s −1 r  s Rcs − Scr Srr Rs Ds Rs .

s =1

Evidently, the two-level Neumann–Neumann preconditioner in Eq. (27) is equipped with a coarse operator R 0T [ F cc ]−1 R 0 to scatter information globally.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

305

4.2.1. Parallel implementation Eq. (11) defining the extended Schur complement system is tackled using PCGM (see Algorithm 1) equipped with the aforementioned two-level Neumann–Neumann preconditioner in Eq. (27). The preconditioned residual vectors at steps 2 and 12 of Algorithm 1 are computed in parallel as outlined in Algorithm 3. Algorithm 3: Two-level Neumann–Neumann preconditioner procedure.

1 2 3

input : (rΓ ) for s = 1, 2, . . . , ns in parallel: do scatter : rΓs = Ds Rs rΓ

compute compute solve update gather

4 5 6 7 8

: : : : :

Frs = Rrs rΓs Fcs = Rcs rΓs Srrs v 1s = Frs s s dcs = Fcs − Scr v1 ns T dc = s=1 Bcs dcs

10

solve : F cc Zc = dc for s = 1, 2, . . . , ns in parallel:

11

do

9

scatter update solve update gather

12 13 14 15 16

output

: : : : :

Zcs = Bcs Zc s v 2s = Frs − Src Zcs s s s Srr Zr = v 2 Z s = Rrs T Zrs + Rcs T Zcs  Z = ns=s 1 RsT Ds Z s

: (Z )

The implementation of Algorithm 3 requires solution of two local fine problems (steps 6 and 14) and a global coarse problem (step 9). The fine problems corresponding to the local stochastic Neumann problems are solved using Algorithm 4. Algorithm 4: Neumann-solver procedure.

1

input

: (Frs )

solve

:

output

: (Zrs )



Aiis Airs s Aris Arr



Xis Zrs



 =

0



Frs

The coarse problem in step 9 of Algorithm 3 has relatively small size. An explicit construction of the coarse Schur complement matrix F cc is computationally intensive. In our implementation, the coarse problem is solved using a parallel PCGM with the Lumped preconditioner as

Mc−1 F cc Zc = Mc−1 dc ,

(28)

where

Mc−1 =

ns

T

s Bcs Acc Bcs .

(29)

s =1

A flow chart of the two-level extended Neumann–Neumann algorithm is shown in Fig. 4. 5. Numerical results The parallel performance of the algorithm is studied for flow through porous media and elasticity problems with spatially varying non-Gaussian material properties. In the PCGM, the forcing term is used as the initial residual and the iterations terminate when

GΓk − SUΓk 2 GΓ0 2

 10−5 .

306

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 4. Flow chart of parallel implementation of the two-level extended Neumann–Neumann algorithm.

Fig. 5. The geometry of the physical domain.

We conducted the numerical investigations in a Linux cluster having 22 nodes (2 Quad-Core 3.0 GHz Intel Xeon processors and 32 GB of memory per node) with InfiniBand interconnect using Message Passing Interface (MPI) [47] and PETSc [46] parallel libraries. The graph partitioning software METIS [58] is used to partition the finite element mesh generated using Gmsh [59]. For visualizations, PMVIS [60], ParaView [61] and MATLAB [62] softwares are used. 5.1. Flow through porous media For numerical illustrations, a two-dimensional steady-state flow through porous media with a random diffusion coefficient is investigated. A lognormal stochastic process is used to model the spatially varying diffusion coefficient, generated from the underlying Gaussian field with an exponential covariance function as

$

2

C αα (x, y) = σ exp

−|x1 − y 1 | b1

+

−|x2 − y 2 | b2

%

.

(30)

For code validation, the results obtained from SSFEM using iterative and direct solvers are compared first (see [43]). Next the accuracy of the results obtained from the proposed parallel iterative algorithm is validated with a serial implementation of the stochastic finite element method as reported in Ref. [63] (see [43]). Moreover, the results obtained from the proposed algorithm are contrasted with the Monte Carlo simulation [43]. The results of the proposed parallel iterative solver for SPDEs and the serial SSFEM code and Monte Carlo simulation are found to be almost identical in the scale of the graph [43]. For numerical investigations, the following values of the parameters are used: σ = 0.3, b1 = 1.0 and b2 = 1.0. The Gaussian random field is assumed to have a constant mean α  = 1.0. The geometry of the physical domain is shown in Fig. 5. This problem is tackled by the authors in [39] using the one-level domain decomposition method. We contrast

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 6. The relative partial sums (

i

j =1

λj/

∞

j =1

307

λ j ) of the sequence of eigenvalues of the exponential covariance function.

Fig. 7. Typical finite element discretization and domain decomposition with 30 subdomains.

the algorithmic performance of the one-level Neumann–Neumann preconditioner (NN) and two-level Neumann–Neumann preconditioner (NNC) to highlight the beneficial effect of the coarse grid. Fig. 6 shows the relative partial sums of the eigenvalues of the exponential covariance function. The first four eigenvalues correspond to about 90% energy of the random field and thus four random variables are used in the KLE of the underlying Gaussian process. The input lognormal diffusion coefficient is represented by four random variables (i.e. four-dimensional) and third order PCE. The corresponding solution process is expanded by the four-dimensional third order PCE. Matthies and Keese [21] (see Theorem 18) proved that the order of PCE for the lognormal diffusion coefficient must be twice the order of PCE for the solution process in order to accurately compute the SSFEM stiffness matrix and the corresponding solution vector. As the primary objective of this investigation is to demonstrate the usefulness of the two-level domain decomposition algorithm, our implementation, using low order PCE approximation of the input lognormal diffusion coefficient, sacrifices some accuracy in favor of computational simplicity in the numerical illustrations. However, the proposed domain decomposition algorithms are general and valid for any arbitrary order of PCE representation for both input and output processes. A typical finite element mesh and its decomposition are sketched in Figs. 7a and 7b respectively. We plot the mean and variance of the solution field in Figs. 8a and 8b respectively. Evidently, these figures highlight the significant influence of uncertainties on the solution field. 5.1.1. Effect of the coefficient of variation For a fixed spatial problem size n = 16,808 and fixed number of subdomains ns = 16, Fig. 9 shows the performance of NN and NNC preconditioners with respect to the magnitude of the coefficient of variation (CoV) of the diffusion coefficient. Although NNC preconditioner requires less number of iterations to converge compared to NN preconditioner, the iteration counts of both the preconditioners remain almost constant as the value of the CoV of the system parameter increases.

308

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 8. The mean and variance of the solution field (adapted from [42]).

Fig. 9. Iteration counts with respect to the magnitude of CoV.

Fig. 10. Iteration counts for different order of PCE.

5.1.2. Effect of the order of PCE For a fixed number of subdomains ns = 128 and increasing spatial problem size, Fig. 10 shows the iteration counts of the NNC preconditioner for different orders of PCE. Evidently, the order of output PCE does not influence the performance of the two-level preconditioner as the number of PCGM iteration is nearly constant.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

309

Fig. 11. The history of the relative PCGM residual.

Fig. 12. Iteration counts versus mesh size for fixed number of subdomains.

5.1.3. Convergence study A typical relative residual decay of PCGM equipped with NN and NNC preconditioners are shown in Fig. 11. Evidently, the PCGM equipped with NNC preconditioner converges in fewer iterations compared to that of the NN preconditioner. 5.1.4. Numerical scalability For the numerical scalability study, we show the PCGM iteration counts while (a) decreasing the finite element mesh size (h) for fixed number of subdomain partitions; (b) increasing number of subdomains (1/ H ) for fixed spatial mesh resolution; (c) increasing the total problem size by including more subdomains with fixed number of elements (H /h). Fig. 12 plots the results for the NN and NNC preconditioners for 128 subdomains and third order PCE while increasing the finite element mesh resolution. Clearly, adding a coarse problem to NN preconditioner significantly reduces the number of iterations for fixed number of subdomains. Increasing problem size does not deteriorate the performance of the NNC preconditioner as the iteration counts remain constant. This fact indicates the numerical scalability of the NNC preconditioner with respect to the problem size. Fig. 13 plots the performance results of the NN and NNC preconditioners for the fixed spatial problem size n = 127,552 using third order PCE with increasing number of subdomains. The NNC preconditioner demands fewer iterations to converge compared to NN preconditioner. The iteration counts for NNC preconditioner remain steady with respect to the number of subdomains. This fact indicates the numerical scalability of the NNC preconditioner in relation to the number of subdomains. For fixed problem size per subdomain n ≈ 8000 and the third order PCE, the addition of new subdomains increases the global problem size as considered next. The corresponding results are shown in Fig. 14. Increasing number of subdomains demonstrates a linear growth in the iteration counts for the NN preconditioner. But the iteration counts remain constant for

310

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 13. Iteration counts versus the number of subdomains for fixed problem size.

Fig. 14. Iteration counts versus the number of subdomains for fixed problem size per subdomain.

the NNC preconditioner. This fact indicates the scalability of the NNC preconditioner in relation to the fixed problem size per subdomain. 5.1.5. Parallel scalability The strong and weak scalability of the NN and NNC preconditioners are plotted in Figs. 15 and 16 respectively. For fixed problem size n = 256,081 (strong scalability), the total execution times for both the NN and NNC preconditioners reduce steadily as expected. For fixed problem size per subdomain n ≈ 8000 (weak scalability) and the third order PCE, the total execution time increases almost linearly in both the preconditioners as we add more processors. The suboptimal weak scaling can be explained as follows. Typically, data should be communicated only among neighboring subdomains as specified by some adjacency matrix. In this case, MPI point-to-point routines (send and receive) [47] should be used to transfer data among the adjacent subdomains. At this stage however, we consider simplicity in our parallel implementation rather than optimality in performance. Hence we handle communication among the subdomains globally via MPI collective routines (scatter and gather) [47]. Although this approach is not an optimal parallel implementation for the proposed domain decomposition algorithm, it can be readily implemented without the need of providing a map of the neighbors to each process and then explicitly managing synchronization among the processes. As the communication among the processes can constitute a substantial part of the execution time, thorough code profiling and tuning (for communication) will be undertaken in the future work since this simple initial parallel implementation looks promising. It is worthwhile to point out that the NNC preconditioner takes longer to execute than the NN preconditioner. This is due to the extra computation incurred in solving the coarse problem at each iteration. Although the total execution time for the PCGM with the NN preconditioner is lower, the PCGM convergence criterion may not be satisfied with the NN preconditioner. The iteration count of the PCGM with the NN preconditioner increases rapidly with the number of subdomains. Furthermore, a scalable preconditioner does not necessarily mean the preconditioner with the minimum computational

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

311

Fig. 15. Strong scalability.

Fig. 16. Weak scalability.

complexity [64], but it ensures convergence while tackling increasingly larger problems. More efficient algorithms that minimize the computational complexities of the coarse problem will reduce the total execution time of the PCGM equipped with the NNC preconditioner. This aspect will be the subject of future investigations. 5.2. Two-dimensional elasticity problem A plane-stress problem is considered next to study the performance of both the NN and NNC preconditioners for twodimensional elasticity problems as shown Fig. 17. For the plane-stress problem, the uncertain constitutive matrix can be represented as (e.g. [12])

⎡

1 ⎣ ν (x, θ) D(x, θ) = 1 − ν 2 (x, θ) 0 E (x, θ)

ν (x, θ) 1 0

0 0

1−ν (x,θ ) 2

⎤

⎦,

(31)

where E (x, θ) and ν (x, θ) are Young’s modulus and Poisson’s ratio respectively. For simplicity and clarity of presentation, we assume Poisson’s ratio to be a spatially invariant deterministic quantity. The spatial uncertainty of the Young’s modulus is modeled as a stochastic process. The methodology for simultaneous treatment of uncertainties in Young’s modules and Poisson’s ratio is described in detail in [65]. However, the proposed algorithm is general and capable of handling arbitrary non-Gaussian models. We represent the Young’s modulus as a lognormal process defined as





E (x, θ) = E 0 (x) exp g (x, θ) ,

(32)

312

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 17. The geometry of the two-dimensional plane-stress problem.

Fig. 18. The relative partial sums (

i

j =1

λj/

∞

j =1

λ j ) of the sequence of eigenvalues of the exponential covariance function.

Fig. 19. Typical finite element discretization and domain decomposition with 30 subdomains.

where E 0 (x) is the mean of Young’s modulus assumed to be 20 × 109 and g (x, θ) is the underlying Gaussian process with a squared exponential covariance function of the following form [32,66]

$ C gg (x, y) = σ 2 exp

−(x2 − x1 )2 − ( y 2 − y 1 )2 b2

% .

(33)

For the numerical investigation, two random variables are used to represent the underlying Gaussian process while the lognormal process is expanded using two-dimensional third order PCE (e.g. [67]). The Gaussian process is assumed to have zero mean and σ and b are 0.35 and 12.65 respectively. In Fig. 18, we show the relative partial sums of the sequence of eigenvalues of the covariance function defined in Eq. (33). 5.2.1. Stochastic characteristics of the solution process The physical domain and its typical decomposition are shown in Figs. 19a and 19b respectively. Figs. 20a and 20b show the mean and variance of the displacement field respectively. The CoV is about 30% at the free end underscoring the effect of uncertainty.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

313

Fig. 20. The mean and variance of the solution process (adapted from [42]).

Fig. 21. Iteration counts versus the order of PCE.

5.2.2. Numerical scalability Fig. 21 demonstrates the performance of both the NN and NNC preconditioners by fixing the problem size and number of subdomains, but increasing the order of PCE of the solution process. Increasing the order of PCE exhibits a linear growth in the iteration counts for the NN preconditioner. However the iteration counts remain steady for the NNC preconditioner. This fact demonstrates that the NNC preconditioner is scalable to the order of PCE. In Fig. 22, the iteration counts for the PCGM with NN and NNC preconditioners, for fixed number of subdomains (ns = 64) and the second order PCE, are plotted while refining the finite element mesh. Evidently, the problem size does not degrade the efficiency of the NNC preconditioner as the number of iterations required for convergence remains the same. The iteration counts increase almost linearly for the NN preconditioner as the problem size grows. The NNC preconditioner takes fewer iterations to converge compared to the NN preconditioner. In Fig. 23, we maintain the problem size in the spatial domain to 111,651 nodes with 223,316 elements but add more subdomains for the second order PCE of the solution process. The results demonstrate that the NN preconditioner is scalable to the number of subdomains. The performance of the NN and NNC preconditioners is shown in Fig. 24 for fixed problem size per subdomain n ≈ 5000. Again these results highlight the fact that the NNC preconditioner is scalable to fixed problem size per subdomain. The performance of the NN preconditioner degrades rapidly as the number of subdomains grows. 5.2.3. Parallel scalability Figs. 25 and 26 demonstrate the strong and weak scalabilities for both NN and NNC preconditioners respectively. For a fixed mesh resolution, both the NN and NNC preconditioners steadily reduce the execution time when we exploit more processors to tackle the problem. For the weak scalability tests, the execution times for both the NN and NNC preconditioners increase linearly. As alluded to previously, the extra computation involved in solving the coarse problem leads to the higher execution time for the PCGM with the NNC preconditioner than that of the NN preconditioner.

314

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

Fig. 22. Iteration counts versus mesh size for fixed number of subdomains.

Fig. 23. Iteration counts versus the number of subdomains for fixed problem size.

Fig. 24. Iteration counts versus the number of subdomains for fixed problem size per subdomain.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

315

Fig. 25. Strong scalability.

Fig. 26. Weak scalability.

6. Conclusions We formulated a novel two-level domain decomposition preconditioner for the iterative substructuring of stochastic PDEs. Within the framework of the intrusive SSFEM, the two-level preconditioner is constructed by adding a coarse grid to the one-level Neumann–Neumann preconditioner. Numerical experiments performed for two-dimensional stochastic problems in flow through porous media and linear elasticity suggest that the probabilistic two-level preconditioner is numerically scalable with respect to the problem size, number of subdomains and fixed problem size per subdomains. The suboptimal weak parallel scalability results can be attributed to the use of collective MPI communication routines to handle the global interface vector. This issue will be addressed in the future investigations by the authors. Acknowledgements The authors gratefully acknowledge the financial support from the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant, Canada Research Chair Program, Canada Foundation for Innovation and Ontario Innovation Trust. References [1] W. Subber, A. Sarkar, Domain decomposition methods of stochastic PDEs: A two-level scalable preconditioner, J. Phys. Conf. Ser. 341 (2012). [2] A. Barth, C. Schwab, N. Zollinger, Multi-level Monte Carlo finite element method for elliptic PDEs with stochastic coefficients, Numer. Math. 119 (2011) 123–161.

316

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

[3] K. Cliffe, M. Giles, R. Scheichl, A. Teckentrup, Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients, Comput. Vis. Sci. 14 (2011) 3–15. [4] D. Xiu, J. Hesthaven, High-order collocation methods for differential equations with random inputs, SIAM J. Sci. Comput. 27 (2005) 1118–1139. [5] I. Babuska, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM J. Numer. Anal. 45 (2007) 1005–1034. [6] F. Nobile, R. Tempone, C. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM J. Numer. Anal. 46 (2008) 2309–2345. [7] J. Back, F.N.L. Tamellini, R. Tempone, Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: A numerical comparison, in: J.S. Hesthaven, E.M. Ronquist (Eds.), Spectral and High Order Methods for Partial Differential Equations, in: Lect. Notes Comput. Sci. Eng., vol. 76, Springer, Berlin, Heidelberg, 2011, pp. 43–62. [8] J. Jakeman, S. Roberts, Stochastic Galerkin and collocation methods for quantifying uncertainty in differential equations: A review, in: G.N. Mercer, A.J. Roberts (Eds.), Proceedings of the 14th Biennial Computational Techniques and Applications Conference, CTAC-2008, ANZIAM J. 50 (2008) C815–C830. [9] H. Elman, C. Miller, E. Phipps, R. Tuminaro, Assessment of collocation and Galerkin approaches to linear diffusion equations with random data, Int. J. Uncertainty Quantif. 1 (2011) 19–33. [10] M. Eldred, J. Burkardt, Comparison of non-intrusive polynomial chaos and stochastic collocation methods for uncertainty quantification, in: 47th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition, 2009. [11] D. Xiu, Efficient collocational approach for parametric uncertainty analysis, Commun. Comput. Phys. 2 (2007) 293–309. [12] R. Ghanem, P. Spanos, Stochastic Finite Element: A Spectral Approach, Springer-Verlag, New York, 1991. [13] R. Ghanem, Ingredients for a general purpose stochastic finite elements implementation, Comput. Methods Appl. Mech. Eng. 168 (1999) 19–34. [14] A. Keese, A review of recent developments in the numerical solution of stochastic partial differential equations (stochastic finite elements), Technical report, Institute of Scientific Computing, Technical University Braunschweig, Brunswick, Germany, 2003, Technical research report: Informatikbericht Nr. 2003-06. [15] M. Schevenels, G. Lombaert, G. Degrande, Application of the stochastic finite element method for Gaussian and non-Gaussian systems, in: Proceeding of ISMA2004, 2004, pp. 3299–3314. [16] O. Kino, O. Maitre, Uncertainty propagation in CFD using polynomial chaos decomposition, Fluid Dyn. Res. 38 (2006) 616–640. [17] R. Ghanem, A. Doostan, On the construction and analysis of stochastic models: Characterization and propagation of the errors associated with limited data, J. Comput. Phys. 217 (2006) 63–81. [18] P. Dostert, Y. Efendiev, T. Hou, Multiscale finite element methods for stochastic porous media flow equations and application to uncertainty quantification, Comput. Methods Appl. Mech. Eng. 197 (2008) 3445–3455. [19] G. Stefanou, The stochastic finite element method: Past, present and future, Comput. Methods Appl. Mech. Eng. 198 (2009) 1031–1051. [20] H. Najm, Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics, Annu. Rev. Fluid Mech. 41 (2009) 35–52. [21] H.G. Matthies, A. Keese, Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations, Comput. Methods Appl. Mech. Eng. 194 (2005) 1295–1331. [22] O. Maitre, O. Knio, Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics, Springer, New York, 2010. [23] I. Babuska, P. Chatzipantelidis, On solving elliptic stochastic partial differential equations, Comput. Methods Appl. Mech. Eng. 191 (2002) 4093–4122. [24] A. Nouy, Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations, Arch. Comput. Methods Eng. 16 (2009) 251–285. [25] D. Xiu, Numerical Methods for Stochastic Computations: A Spectral Method Approach, Princeton University Press, New Jersey, 2010. [26] E. Rosseel, S. Vandewalle, Iterative solvers for the stochastic finite element method, SIAM J. Sci. Comput. 32 (2010) 372–397. [27] R. Ghanem, R. Kruger, Numerical solution of spectral stochastic finite element systems, Comput. Methods Appl. Mech. Eng. 129 (1996) 289–303. [28] M. Pellissetti, R. Ghanem, Iterative solution of systems of linear equations arising in the context of stochastic finite elements, Adv. Eng. Softw. 31 (2000) 607–616. [29] A. Keese, H. Matthies, Parallel computation of stochastic groundwater flow, in: NIC Symposium, vol. 20, John von Neumann Institute for Computing, 2004, pp. 399–408. [30] C. Powell, H. Elman, Block-diagonal preconditioning for spectral stochastic finite-element systems, IMA J. Numer. Anal. 29 (2009) 350–375. [31] D. Ghosh, P. Avery, C. Farhat, A FETI-preconditioned conjugate gradient method for large-scale stochastic finite element problems, Int. J. Numer. Methods Eng. 80 (2009) 914–931. [32] E. Ullmann, A Kronecker product preconditioner for stochastic Galerkin finite element discretizations, SIAM J. Sci. Comput. 32 (2010) 923–946. [33] O. Maitre, O. Knio, B. Debusschere, H. Najm, R. Ghanem, A multigrid solver for two-dimensional stochastic diffusion equations, Comput. Methods Appl. Mech. Eng. 192 (2003) 4723–4744. [34] H. Elman, D. Furnival, Solving the stochastic steady-state diffusion problem using multigrid, IMA J. Numer. Anal. 27 (2007) 675–688. [35] E. Rosseel, T. Boonen, S. Vandewalle, Algebraic multigrid for stationary and time-dependent partial differential equations with stochastic coefficients, Numer. Linear Algebra Appl. 15 (2008) 141–163. [36] D. Chung, M. Gutiérrez, L. Graham-Brady, F. Lingen, Efficient numerical strategies for spectral stochastic finite element models, Int. J. Numer. Methods Eng. 64 (2005) 1334–1349. [37] A. Sarkar, N. Benabbou, R. Ghanem, Domain decomposition of stochastic PDEs: Theoretical formulations, Int. J. Numer. Methods Eng. 77 (2009) 689–701. [38] W. Subber, A. Sarkar, Domain decomposition of stochastic PDEs: A novel preconditioner and its parallel performance, in: High Performance Computing Symposium, in: Lect. Notes Comput. Sci., vol. 5976, Springer, 2010, pp. 251–268. [39] W. Subber, A. Sarkar, Performances of one-level and two-level primal domain decomposition preconditioners for stochastic PDEs, Probab. Eng. Mech. (2013), submitted for publication. [40] W. Subber, A. Sarkar, Primal and dual–primal iterative substructuring methods of stochastic PDEs, J. Phys. Conf. Ser. 256 (2010). [41] W. Subber, A. Sarkar, Domain decomposition methods of stochastic PDEs, in: R. Bank, M. Holst, O. Widlund, J. Xu (Eds.), Domain Decomposition Methods in Science and Engineering XX, in: Lect. Notes Comput. Sci. Eng., vol. 91, Springer, Berlin, Heidelberg, 2013, pp. 191–198. [42] W. Subber, A. Sarkar, Dual–primal domain decomposition method for uncertainty quantification, Comput. Methods Appl. Mech. Eng. 266 (2013) 112–124. [43] W. Subber, Domain decomposition methods for uncertainty quantification, PhD thesis, Department of Civil and Environmental Engineering, Carleton University, 2012. [44] J. Cros, A preconditioner for the Schur complement domain decomposition method, in: Fourteenth International Conference on Domain Decomposition Methods, National Autonomous University of Mexico (UNAM), 2003, pp. 373–380. [45] B. Sousedik, J. Mandel, On the equivalence of primal and dual substructuring preconditioners, Electron. Trans. Numer. Anal. 31 (2008) 384–402. [46] S. Balay, K. Buschelman, W. Gropp, D. Kaushik, M. Knepley, L. McInnes, B. Smith, H. Zhang, PETSc web page http://www.mcs.anl.gov/petsc, 2009. [47] Message passing interface forum, http://www.mpi-forum.org, 2009. [48] W. Subber, H. Monajemi, M. Khalil, A. Sarkar, A scalable parallel uncertainty analysis and data assimilation framework applied to some hydraulic problems, in: International Symposium on Uncertainties in Hydrologic and Hydraulic, Montreal, Canada, 2008.

W. Subber, A. Sarkar / Journal of Computational Physics 257 (2014) 298–317

317

[49] A. Toselli, O. Widlund, Domain Decomposition Methods – Algorithms and Theory, Springer Ser. Comput. Math., vol. 34, Springer, Berlin, 2005. [50] A. Quarteroni, A. Valli, Domain Decomposition Methods for Partial Differential Equations, Numerical Mathematics and Scientific Computation, Oxford University Press, USA, 1999. [51] B. Smith, P. Bjorstad, W. Gropp, Domain Decomposition: Parallel Multilevel Methods for Elliptic Partial Differential Equations, Cambridge University Press, New York, 1996. [52] T. Mathew, Domain Decomposition Methods for the Numerical Solution of Partial Differential Equations, Lect. Notes Comput. Sci. Eng., vol. 61, Springer, Berlin, 2008. [53] Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM, Philadelphia, 2003. [54] C. Kelley, Iterative Methods for Linear and Nonlinear Equations, SIAM, Philadelphia, 1995. [55] M. Benzi, Preconditioning techniques for large linear systems: A survey, J. Comput. Phys. 182 (2002) 418–477. [56] C. Farhat, M. Lesoinne, K. Pierson, A scalable dual–primal domain decomposition method, Numer. Linear Algebra Appl. 7 (2000) 687–714. [57] J. Mandel, C. Dohrmann, R. Tezaur, An algebraic theory for primal and dual substructuring methods by constraints, Appl. Numer. Math. 54 (2005) 167–193. [58] G. Karypis, V. Kumar, METIS unstructured graph partitioning and sparse matrix ordering system, University of Minnesota, Department of Computer Science and Engineering, 1995, http://glaros.dtc.umn.edu/gkhome/metis/metis/overview. [59] J.R.C. Geuzaine, Gmsh: A three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Methods Eng. 79 (2009) 1309–1331. [60] B. Oztekin, PMVIS partitioned mesh visualizer, University of Minnesota, Department of Computer Science and Engineering, 2012, http://www-users.cs. umn.edu/~oztekin/pmvis/. [61] ParaView, http://www.paraview.org, 2012. [62] MATLAB, http://www.mathworks.com/products/matlab, 2012. [63] P. Frauenfelder, C. Schwab, R. Todor, Finite elements for elliptic problems with stochastic coefficients, Comput. Methods Appl. Mech. Eng. 194 (2005) 205–228. [64] T. Chan, T. Mathew, Domain decomposition algorithms, Acta Numer. 3 (1994) 61–143. [65] R. Ghanem, G. Saad, A. Doostan, Efficient solution of stochastic systems: Application to the embankment dam problem, Struct. Saf. 29 (2007) 238–251, A benchmark study on reliability in high dimensions. [66] D. Griffiths, G. Fenton (Eds.), Probabilistic Methods in Geotechnical Engineering, CISM Courses and Lectures, vol. 491, Springer, 2007. [67] R. Ghanem, The nonlinear Gaussian spectrum of log-normal stochastic processes and variables, J. Appl. Mech. 66 (1999) 964–973.