Basis adaptation and domain decomposition for steady-state partial differential equations with random coefficients

Journal of Computational Physics 351 (2017) 203–215

Contents lists available at ScienceDirect

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

Basis adaptation and domain decomposition for steady-state partial differential equations with random coefficients R. Tipireddy, P. Stinis, A.M. Tartakovsky ∗ Pacific Northwest National Laboratory, Richland, WA, USA

a r t i c l e

i n f o

Article history: Received 21 December 2015 Received in revised form 18 February 2017 Accepted 31 August 2017 Available online 4 September 2017 Keywords: Basis adaptation Dimension reduction Domain decomposition Polynomial chaos Uncertainty quantification

a b s t r a c t We present a novel approach for solving steady-state stochastic partial differential equations in high-dimensional random parameter space. The proposed approach combines spatial domain decomposition with basis adaptation for each subdomain. The basis adaptation is used to address the curse of dimensionality by constructing an accurate low-dimensional representation of the stochastic PDE solution (probability density function and/or its leading statistical moments) in each subdomain. Restricting the basis adaptation to a specific subdomain affords finding a locally accurate solution. Then, the solutions from all of the subdomains are stitched together to provide a global solution. We support our construction with numerical experiments for a steady-state diffusion equation with a random spatially dependent coefficient. Our results show that accurate global solutions can be obtained with significantly reduced computational costs. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Uncertainty quantification for systems with a large number of deterministically unknown input parameters is a formidable computational task. In the probabilistic framework, uncertain parameters are treated as random variables/fields, yielding the governing equations stochastic. The two most popular methods for solving stochastic equations are Monte Carlo (MC) and polynomial chaos (PC) expansion. In both methods, the random input parameters are represented with d random variables using truncated Karhunen–Loève (KL) expansion [1]. MC methods are robust and easy to implement, but they converge at a very slow rate. Hence, they require a large number of samples. On the other hand, the MC convergence rate does not depend on d. Contrary to this, the computational cost of standard PC methods increases exponentially with increasing d, a phenomenon often referred to as the “curse of dimensionality” [2–5]. Because of this, standard PC methods are only efficient for small to moderate d [6–8]. Recently, basis adaptation [9] and active subspace methods [10] have been presented to identify a low-dimensional representation of the solution of stochastic equations. In the current work, we propose a new approach that combines basis adaptation with spatial domain decomposition ([11–15]) to address the problem posed by large d. In particular, we focus on partial differential equations (PDEs) with spatially dependent random coefficients. If considered in the whole spatial domain, the randomness is very high-dimensional. However, a few dominant random parameters could suffice locally [16]. As a result, we can obtain low-dimensional local representations of the random solution in each subdomain. To obtain the low dimensional representation in each subdomain, we use the Hilbert space KL expansion [17]. Then, we reconstruct a global solution by stitching together the local solutions for each subdomain.

*

Corresponding author. E-mail address: [email protected] (A.M. Tartakovsky).

http://dx.doi.org/10.1016/j.jcp.2017.08.067 0021-9991/© 2017 Elsevier Inc. All rights reserved.

204

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

The paper is organized as follows: Section 2 presents the problem of uncertainty quantification for PDEs with random coefficients. Section 3 introduces our approach that combines the domain decomposition and basis adaptation. Section 4 contains numerical results for a two-dimensional steady-state diffusion equation for two different types of boundary conditions. Section 5 presents conclusions and ideas for future work. 2. PDEs with random coefficients Let D be an open subset of Rn and (, , P ) be a complete probability space with sample space , probability measure P . We want to find a random field, u (x, ω) : D ×  → R, such that P -almost surely

σ -algebra  , and

L(x, u (x, ω); a(x, ω)) = f (x, ω) in D × , B(x, u (x, ω); a(x, ω)) = h(x, ω) on ∂ D × ,

(1)

where L is a differential operator and B is a boundary operator. We model the uncertainty in the stochastic PDE (1) by treating the coefficient a(x, ω) as a random field and compute the effect of this uncertainty on the solution field u (x, ω). To solve the stochastic PDE numerically, we discretize the random fields a(x, ω) and u (x, ω) in both spatial and stochastic domains. In this paper, we will focus on the special case where a(x, ω) is a log-normal random field [18]. This means that a(x, ω) = exp[ g (x, ω)], where g (x, ω) is a Gaussian random field whose mean and covariance function are known. We approximate g (x, ω) with a truncated KL expansion, while the coefficient a(x, ω) and solution u (x, ω) are approximated through truncated PC expansions. 2.1. Karhunen–Loève expansion of the random field g (x, ω) The random field g (x, ω) can be approximated with a truncated KL expansion [1],

g (x, ω) ≈ g (x, ξ (ω)) = g 0 (x) +

d  

λi g i (x)ξi (ω),

(2)

i =1

where d is the number of random variables in the truncated expansion; ξ = (ξ1 , · · · , ξd ) T , ξi are uncorrelated random variables with zero mean; g 0 (x) is the mean of the random field g (x, ω); and (λi , g i (x)) are eigenvalues and eigenvectors obtained by solving the eigenvalue problem,

 C g (x1 , x2 ) g i (x2 )dx2 = λi g i (x1 ),

(3)

D

where C g (x1 , x2 ) is the covariance function of the Gaussian random field g (x, ω). The eigenvalues are positive and nonincreasing, and the eigenfunctions g i (x) are orthonormal,



g i (x) g j (x)dx = δi j ,

(4)

D

where δi j is the Kronecker delta. In this work, we assume the random variables ξi have Gaussian distribution. Therefore, ξi are independent. 2.2. Polynomial chaos expansion We approximate the input random field, a(x, ω) and the solution field u (x, ω) using truncated PC expansions [19] in Gaussian random variables as follows:

a(x, ω) ≈ a(x, ξ (ω)) = a0 (x) +

Nξ 

ai (x)ψi (ξ )

(5)

i =1

and

u (x, ω) ≈ u (x, ξ (ω)) ≈ u 0 (x) +

Nξ 

u i (x)ψi (ξ ),

(6)

i =1

(d+ p )!

where N ξ = d! p ! − 1 is the number of terms in PC expansion for dimension d and order p, u 0 (x) is the mean of the solution field, u i (x) are PC coefficients, and {ψi (ξ )} are multivariate Hermite polynomials. These polynomials are orthogonal with respect to the inner product defined by the expectation in the stochastic space,

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

205

 ψi (ξ ), ψ j (ξ ) ≡

ψi (ξ (ω))ψ j (ξ (ω))d P (ω) = δi j ,

(7)



where d P (ω) is the diagonal d-dimensional Gaussian distribution. Once the random coefficient a(x, ω) and the solution field u (x, ω) are approximated, the stochastic PDE (1) is transformed into

L(x, ξ , u (x, ξ ); a(x, ξ )) = f (x, ξ ) in D × , B(x, ξ , u (x, ξ ); a(x, ξ )) = h(x, ξ ) on ∂ D × .

(8)

Intrusive methods, such as stochastic Galerkin [2], or non-intrusive methods, such, as sparse-grid collocation [3], can be used to solve the parameterized stochastic PDE (8). However, both methods suffer from the curse of dimensionality, i.e., their computational costs exponentially increase with increasing d and/or p. To address the curse of dimensionality, we propose a novel approach based on the basis adaptation method, originally introduced in [9] to compute locally accurate solutions. Our approach combines the basis adaptation method with spatial domain decomposition, which allows us to compute an accurate solution in the entire domain. In the following section, we demonstrate how basis adaptation can be applied to find a low-dimensional solution in non-overlapping spatial subdomains. 3. Basis adaptation in a spatial subdomain To address the curse of dimensionality in solving (1), we first decompose the spatial domain D into a set of nonoverlapping subdomains. Then, we use the basis adaptation to find a low-dimensional solution space in each subdomain. 3.1. Domain decomposition We divide the spatial domain D ⊂ Rn into a set of non-overlapping subdomains [16], D s ⊂ D , s = 1, · · · , S, such that

D=

S 

Ds,

D s ∩ D s = ∅, s = s .

(9)

s =1 s In each subdomain D s , we find a low-dimensional stochastic basis (PC) for a new set of random variables η˜ = {η1s , · · · , ηrs }, r d, adapted to that subdomain, such that the solution can be well represented in the low-dimensional

basis. When the stochastic basis is adapted to the subdomain D s , the low-dimensional solution computed by solving (1) is accurate only in D s . Hence, we repeat the process in each subdomain to get an accurate low-dimensional solution in the entire spatial domain. 3.2. Basis adaptation through Hilbert–Karhunen–Loève expansion

Here, we extend the basis adaptation method of [9] to construct a low-dimensional solution representation in each subdomain D s using the Hilbert–Karhunen–Loève expansion. The KL expansion [1] described in Section 2.1 provides optimal representation of a random field in L 2 space. However, to take advantage of this property for PDEs with random coefficients, the solution should satisfy certain regularity and smoothness conditions [17]. Thus, it is appropriate to seek a representation of the solution field u (x, ω) in a subset of L 2 (). Such a restricted expansion has been explored in various applications [20–24]. In this work, we first compute the Gaussian part (up to linear terms in PC expansion) of the solution in the whole domain D, i.e.,

u g (x, ξ (ω)) = u 0 (x) +

d 

u i (x)ξi .

(10)

i =1

Note that in realistic situations, this computation could be feasible even if the full solution computation is not. Let u sg (x, ξ (ω)) ⊂ u g (x, ξ (ω)) be the Gaussian part of the solution in subdomain D s , i.e.,

u sg (x, ξ (ω)) = u g (x, ξ (ω))I D s (x),

(11)

where I is the indicator function such that for any set α (Iα = 1, if x ∈ α , and Iα = 0, if x ∈ / α ). For each subdomain D s , we construct the covariance function of the Gaussian part of the solution u sg (x, ξ ) as follows:

C us g (x1 , x2 ) =

d 

u i (x1 )u i (x2 ),

x1 , x2 ∈ D s .

i =1

The Hilbert space KL expansion of u sg (x, ξ ) (see e.g., [17]) is

(12)

206

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

u sg (x, ξ (ω)) = u 0s (x) +

d  

μis φis (x)ηis (ω), x ∈ D s

(13)

i =1

where (μis , φis (x)) are eigenpairs corresponding to the above expansion. They can be obtained by solving the eigenvalue problem:

 C us g (x1 , x2 )φis (x1 )dx1 = μis φis (x2 ),

i = 1, 2, · · · .

(14)

Ds

After solving the eigenvalue problem, we can use (13) to write the random variable



1

ηis = 

μ

(u sg (x, ξ ) − u˜ 0s (x))φis (x)dx,

s i Ds



1

=

μis D

ηis as

(u 0 (x) +

d 

u j (x)ξ j − u˜ 0s (x))φis (x)dx.

i = 1, 2, · · · , d

(15)

j =1

s

Note that u 0 (x) = u 0s (x) by construction. Hence, s i

η =

⎛ ⎞ d  ⎝ u j (x)ξ j ⎠ φ s (x)dx,



1

μis D ⎛

i

x ∈ D s , i = 1, 2, · · · , d

j =1

s

⎞  d  1 ⎜ ⎟ = u j (x)φis (x)dx⎠ ξ j , ⎝

μis D

j =1

x ∈ D s , i = 1, 2, · · · , d

s

= aisj ξ j , where aisj = 1 s μi

(16)

Ds

u j (x)φis (x)dx, i , j = 1, · · · , d. This provides a linear relation between

ηis and {ξ j }. Because, ξ j are stan-

dard Gaussian random variables, the new variables {ηis } also are Gaussian random variables and can be normalized to obtain standard Gaussian random variables. The resulting Gaussian random variables {ηis } are uncorrelated and, thus, also independent. We can reformulate the stochastic PDE (1) in terms of {ηis } and solve it using intrusive or non-intrusive methods. 3.3. Dimension reduction in the subdomain D s ⊂ D Let A s = [aisj ] be an isometry in Rd and define

ηs = A s ξ ,

η s as

As As T = I ,

s 1, · · ·

(17)

where η = {η ,η is a vector of standard normal random variables. With the mapping in (17), Gaussian Hilbert space spanned by random variables ξ in D s . This means s

s T } d

N

u (x, ξ (ω)) = u

As

s

(x, η (ω)) =

A u 0 s (x) +

s

η is defined in the

s

η 

A

u i s (x)ψi (η s ),

x ∈ Ds.

(18)

i =1

For this construction to be worth the effort, we should be able to obtain an accurate representation of the solution without having to use functions of all d new variables ηis , i = 1, . . . , d. This depends on how fast the magnitude of the eigenvalues μis (i = 1, . . . , d) decays. For the steady-state diffusion equation considered here, the solution is expected to be smoother than the input random field a(x, ω). In addition, we employ the Hilbert space KL expansion in a spatial subdomain D s ⊂ D. Thus, eigenvalues μis in subdomain D s are expected to decay faster than the eigenvalues λi of the input random field (see Fig. 2). Because of this, it is sufficient to perform the Hilbert space KL expansion in (13) only in r( d) dimensions. s With the new set of random variables η˜ = {η1 , · · · , ηr }, (8) can be reformulated as

L(x, ξ , η˜ s , u (x, η˜ s ); a(x, ξ )) = f (x, ξ ) in D × , B(x, ξ , η˜ s , u (x, η˜ s ); a(x, ξ )) = h(x, ξ ) on ∂ D × .

(19)

The stochastic PDE (19) can be solved using non-intrusive or intrusive methods (described in the two sections that follow).

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

207

3.4. Sub-domain solution using intrusive methods s

Expanding the solution of (19) with Hermite polynomials in η˜ and the random coefficient a(x, ξ ) with Hermite polynomials in ξ , and performing a Galerkin projection leads to the deterministic PDEs, N

η Nξ   s

A

s

s

Li [u˜ j s ] E [ψi (ξ )ψ j (η˜ )ψk (η˜ )] = 0,

k = 1, · · · , N η s ,

(20)

j =0 i =0 A

for the expansion coefficients u˜ j s , j = 1, . . . , N η s . Also, in (20), Li = (L, ψi (ξ )) is the inner product of the operator L with

the PC coefficient ψi (ξ ). If the number of retained random variables ηis is small, then the number of PC terms N η s is much smaller than N ξ . Although we find a low-dimensional stochastic basis adapted to a subdomain D s , we solve (20) in the entire spatial domain D. However, because the random variables ηis are adapted to the subdomain D s , the solution is accurate only in that subdomain. Hence, for each subdomain D s , s ∈ S, we solve (20) with stochastic basis adapted to that subdomain and retain the solution corresponding to that subdomain while discarding the solution for the remaining subdomains. In this way, we obtain a solution for each subdomain using the adaptive basis in that subdomain. Let | S | be the total number of non-overlapping subdomains. If the subdomains are of equal size (and resolution), the computational cost for obtaining the solution in the entire domain is | S | times the cost of solving (20) plus the cost of solving (8) up to the first order PC. 3.5. Sub-domain solution using non-intrusive methods In a non-intrusive method, a PDE is solved at predefined quadrature points (in random space). We use sparse-grid collocation points based on Smolyak approximation [25,26]. Unlike the tensor–product of one dimensional quadrature points, the sparse-grid method judiciously chooses products with only a small number of quadrature points. These product rules depend on an integer value called sparse-grid level (see [25,26]). As the order of PC expansion increases, we need higher sparse-grid levels to maintain solution accuracy. In the sparse-grid method, the number of collocation points increases with s the stochastic dimension (number of random variables) and the sparse-grid level. Let η˜ q be the collocation points associated s with the random variables η˜ . Then, the corresponding points ξ qs are obtained as

1 ˜ qs , ξ qs = [ A − s ]r η

(21)

1 −1 where the d × r matrix [ A − s ]r consists of the first r columns of A s . Finally, the deterministic PDE

L[u (ξ qs )] = 0

(22) As

is solved for each collocation point ξ qs (q = 1, · · · , Q ), and the PC coefficient u˜ i A

u˜ i s (x) =

Q 

is computed by projection,

s

u (ξ qs )ψi (η˜ q ) w q ,

(23)

q =1

where w q are weights for the quadrature points. s

3.6. The error of the low-dimensional solution, u˜ A s (x, η˜ ) In the basis adaptation method, if the dimension of the new basis in η s and the dimension of original basis in ξ are equal, there is no dimension reduction, and the accuracy of the solution u˜ A s (x, η s ) in terms of η s is the same as that of the solution u (x, ξ ). This means:

 (x, ω) = u (x, ξ (ω)) − u˜ A s (x, η s (ω)) = 0, if r = d.

(24)

However, if the eigenvalues μis decay fast (see Fig. 2), we can truncate the Hilbert space KL expansion in (13) to r + 1 terms and represent the solution in a low-dimensional (r d) space, A

u (x, ξ (ω)) = u˜ 0 s (x) +



s

u˜ I1 ψI1 (η˜ ) +

I1

⎛

u (x, ξ (ω)) − ⎝u˜ 0 s (x) + A

 I1



s

s

u˜ I12 ψI12 (η˜ , ηˆ ) +

I12

⎞

u˜ I1 ψI1 (η˜ )⎠ = s



s

u˜ I2 ψI2 (ηˆ ),

I2

 I12

s

s

u˜ I12 ψI12 (η˜ , ηˆ ) +

 I2

s

u˜ I2 ψI2 (ηˆ ),

208

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

s

u (x, ξ (ω)) − u˜ A s (x, η˜ (ω)) =

 (x, ω) =



 I12

s

s

s

u˜ I12 ψI12 (η˜ , ηˆ ) +

s

u˜ I12 ψI12 (η˜ , ηˆ ) +

I12





s

u˜ I2 ψI2 (ηˆ ),

I2 s

u˜ I2 ψI2 (ηˆ ),

if r d

(25)

I2 s

s

s s where η˜ = {η1s , · · · , ηrs }, ηˆ = {ηrs+1 , · · · , ηds }, multi-indices I1 and I2 correspond to PC expansion terms in η˜ and ηˆ , s s respectively, and I12 is the multi-index for mixed PC terms in η˜ and ηˆ .

4. Numerical examples In this section, we apply the proposed approach to a two-dimensional steady-state diffusion equation with random coefficient, defined on the spatial domain D = [0, 240] × [0, 60] such that x = (x1 , x2 ) ∈ D. We solve this equation for two types of boundary conditions. In the first case, we use Dirichlet boundary conditions at the boundaries perpendicular to the x1 direction and Neumann boundary conditions at the other two boundaries. In the second case, we use Dirichlet boundary conditions on all four sides. Let the random coefficient a(x, ω) : D ×  → R be bounded and strictly positive,

0 < al ≤ a(x, ω) ≤ au < ∞ a.e.

in

D × .

(26)

In the first case, we solve the boundary-value problem:

−∇ · (a(x, ω)∇ u (x, ω)) = f (x, ω) in D × , u (x, ω) = 100 on x1 = 0, u (x, ω) = 10 on x1 = 240,

 · ∇ u (x, ω) = 0 on x2 = 0, n  · ∇ u (x, ω) = 0 on x2 = 60, n

(27)

where u (x, ω) : D ×  → R. In the second case, we solve the boundary-value problem:

−∇ · (a(x, ω)∇ u (x, ω)) = f (x, ω) in D × , u (x, ω) = 100 on x1 = 0, u (x, ω) = 10 on x1 = 240, u (x, ω) = 0 on x2 = 0, u (x, ω) = 0 on x2 = 60.

(28)

We assume that random coefficient a(x, ω) = exp[ g (x, ω)] has log-normal distribution with mean a0 (x) = 5.0 and standard deviation σa = 2.5. The Gaussian random field g (x, ω) has correlation function



C g (x, y ) = σ

2 g exp

−

(x1 − y 1 )2 l21

−

 

where the standard deviation

(x2 − y 2 )2

σ g = ln 1 + a

l22 σa

2 0 (x)





in

D × ,

(29)

 , the mean g 0 (x) = ln

  a0 (x) ( 1+ σa 2 ) a0 (x)

, and the correlation lengths l1 = 24

and l2 = 20. We decompose the spatial domain into eight subdomains and independently find a low-dimensional solution in each subdomain using the KL expansion and basis adaptation methods described in Sections 3.2 and 3.3. The number of stochastic dimensions d of the problem (1) and (2) depends on the decay rate of eigenvalues in the g (x, ω) KL expansion, which are shown in Fig. 2. Here, for computational considerations, we assume that the first 10 eigenvalues are sufficient to capture the uncertainty in a(x, ω) = exp[ g (x, ω)] and use 10 dimensions to describe the reference system. To evaluate the accuracy of the solution obtained with the reduced model, we solve the full stochastic system in the domain D with stochastic dimension 10 and sparse-grid level 5. This requires Q = 8761 collocation points. In the domain decomposition with basis adaptation approach, we initially solve the full system in the entire domain D with a coarse grid (sparse-grid level 3) and dimension 10, which required Q = 221 collocation points. We use a coarser grid to obtain the Gaussian component of the solution and identify the low-dimensional space directions in each spatial subdomain. To find the reduced dimension in each subdomain, we use KL expansion (13) of the Gaussian solution u sg (x, ω). We truncate this KL expansion after 3 eigenvalues because starting with the third eigenvalue, the eigenvalues of u sg (x, ω) are more than one order of magnitude smaller than the corresponding eigenvalues of g (x, ω) (see Fig. 2). Fig. 2 also shows that the eigenvalues for different subdomains decay uniformly, hence we use a dimension of 3 for each subdomain. The low-dimensional representation requires 165 collocation points (dimension = 3 and sparse-grid level 5) per subdomain. As result, the total number the

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

209

Fig. 1. Spatial domain decomposed into eight non-overlapping subdomains.

Fig. 2. Decay of the eigenvalues of covariance function of g (x, ω) in domain D and the covariance function of Gaussian solution in subdomains D 1 , D 2 , D 3 , and D 4 .

Fig. 3. Mean of u in subdomain D 1 : (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

corresponding deterministic PDE (22) that must be solved in the domain decomposition with basis adaptation method is 8 × 165 + 221 = 1541, which is much smaller than Q = 8761 needed to solve the problem without basis adaptation. Fig. 1 depicts the decomposition of the spatial domain into eight subdomains. Fig. 2 compares the decay of eigenvalues of the covariance function of g (x, ω) in domain D and those of the covariance function of the Gaussian solution u g (x, ω) in subdomains D 1 , D 2 , D 3 , and D 4 . This figure demonstrates that eigenvalues decay fast in all subdomains, which justifies using a low-dimensional solution in each subdomain after basis adaptation. Figs. 3–6 present the mean and standard deviation in subdomains D 1 and D 2 , computed with full dimension (d = 10) and reduced dimension (r = 3) for the first type of boundary conditions. For example, Fig. 3a corresponds to the mean of the solution computed with full dimension d = 10 and sparse-grid level 5. Fig. 3b shows the mean of the solution computed with the low-dimensional r = 3 basis adapted to the subdomain D 1 . Fig. 3c shows the error in the mean computed as the difference between quantities in Fig. 3a and Fig. 3b. Note that although the basis is adapted only to subdomain D 1 , the low-dimensional solution is computed in the entire spatial domain D. Hence, in Fig. 3, we are only interested in the solution in subdomain D 1 . We can see the error in subdomain D 1 is quite small compared to other domains. Similarly, Fig. 4 depicts a comparison of the mean solution in subdomain D 2 .

210

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

Fig. 4. Mean of u in subdomain D 2 : (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 5. Standard deviation of u in subdomain D 1 : (a) Standard deviation of u (x, ξ ) with d = 10 and p = 3; (b) standard deviation of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 6. Standard deviation of u in subdomain D 2 : (a) Standard deviation of u (x, ξ ) with d = 10 and p = 3; (b) standard deviation of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

A similar comparison for the standard deviation is shown in Figs. 5 and 6. We also computed and collected the adapted solution from each subdomain and plotted the mean and standard deviations in Figs. 7 and 8, respectively. For quantitative comparison, we compute the relative error in mean as,

μe =

μu (x) − μu˜ (x)2 , μu (x)2

(30)

where μu (x) is the mean of the full solution u (x, ξ ) and μu˜ (x) is the mean of the reduced solution u (x, η ) computed using the basis adaption and domain decomposition method and  · 2 is the L 2 norm. Similarly, the relative error in standard deviation is

σe =

σu (x) − σu˜ (x)2 , σu (x)2

(31)

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

211

Fig. 7. Mean of u in entire domain D for the boundary conditions in (27): (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 8. Standard deviation of u in entire domain D for the boundary conditions in (27): (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 9. Mean of u in entire domain D for the boundary conditions in (28): (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

where σu (x) and σu˜ (x) are the standard deviation of u (x, ξ ) and u (x, η ), respectively. Relative errors in the mean and standard deviation in Figs. 7 and 8 are μe = 0.0087 and σe = 0.0281, respectively, which demonstrates robustness of the proposed method. For the second case (Dirichlet boundary conditions on all sides of the domain), we show the mean and standard deviation in Figs. 9 and 10. Relative errors in mean and standard deviation in these figures are μe = 0.0004 and σe = 0.0078, respectively. We observe very good agreement for the mean and standard deviation results obtained with full and reduced dimensions. We also recognize a small mismatch of the solutions from adjacent subdomains at the subdomain boundaries. This discrepancy is due to the low-dimensional approximation, which uses a different adapted basis in each subdomain. If required, this inconsistency can be eliminated by appropriate post-processing, such as averaging, interpolation, or projecting each low-dimensional solution onto a global basis.

212

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

Fig. 10. Standard deviation of u in entire domain D for the boundary conditions in (28): (a) Mean of u (x, ξ ) with d = 10 and p = 3; (b) mean of u (x, η) with r = 3 and p = 3; and (c) absolute error between results in (a) and (b). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 11. Probability density function of u (x, ξ ) and u (x, η), obtained with d = 10, r = 3, and p = 3, at points: a) (24, 15) in D 1 ; b) (81, 15) in D 2 ; and c) (150, 15) in D 3 .

Fig. 12. Probability density function of u (x, ξ ) and u (x, η), obtained with d = 10, r = 3, and p = 3, at points: a) (210, 15) in D 4 ; b) (210, 45) in D 5 ; and c) (150, 45) in D 6 .

A complete solution of a stochastic PDE is given by the probability density function (PDF). Figs. 11, 12, and 13 show PDF of u (x, ξ ) and u (x, η ) computed in each subdomain using the full stochastic dimension and the reduced dimension adapted to that subdomain, respectively. A good agreement is observed between high- and low-dimensional PDF estimates. Although the local solution u (x, η s ) obtained from basis adaptation in a subdomain D s is expected to be accurate in

only that subdomain, cross correlation of the solution across adjacent subdomains (between u (x, η s ) and u (x, η s )) is still s preserved. This is because the local basis η in each subdomain is a linear function of the global basis ξ and hence the local basis across adjacent domains is correlated. To demonstrate this, we compute the covariance function of u in full dimension and in reduced dimension from local solutions. Fig. 14a shows the covariance function computed from the full solution u (x, ξ ) in terms of global basis and Fig. 14b shows that computed from the reduced solution u (x, η s ) in terms of local basis in each subdomain. Fig. 14c shows the error as a difference between the two quantities. We can see that the global covariance function is recovered from local solutions. We also compare the eigenvalues and eigenfunctions of the covariance functions of u (x, ξ ) and u (x, η ). Fig. 15 shows a very good agreement between eigenvalues of these covariance functions. In Fig. 16, we compare first three eigenfunctions of the covariance functions. These results show that the low dimensional solution computed in each subdomain is an accurate representation of the global solution and, when all local solutions are stitched together, it accurately represents the global solution in entire domain with correct statistical properties.

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

213

Fig. 13. Probability density function of u (x, ξ ) and u (x, η ), obtained with d = 10, r = 3, and p = 3, at points: a) (81, 45) in D 7 ; and b) (24, 45) in D 8 .

Fig. 14. Covariance function of: a) u (x, ξ ), d = 10, b) u (x, η), r = 3, c) error between a) and b). Results are obtained in the entire domain D for the boundary conditions in (27). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Fig. 15. Eigenvalues of the covariance function of u (x, ξ ), d = 10 (in black) and u (x, η ), r = 3 (in red) for the boundary conditions in (27).

5. Conclusions We have presented a novel approach to accelerate solution of time-independent PDEs with random parameters. The main idea of the proposed approach is to decompose the spatial domain into a set of non-overlapping subdomains and use a stochastic basis adaptation method to compute a low-dimensional representation of the solution in each subdomain. We have provided numerical results where we compared the mean, standard deviation, and PDF of the state variable u computed both with full and reduced dimensional representations. In our numerical experiments, the low-dimensional solutions agree quite well with high-dimensional solutions. Moreover, they do so at a significantly smaller computational cost. In the current work, we employ the simplest method to obtain a global solution by stitching together the solutions from the subdomains. As expected, this creates spikes in the error along subdomain boundaries, because the solution in each subdomain is computed separately in a different low dimensional space. For the cases examined here, these spikes in the error are relatively small. However, they can be eliminated by using an appropriate interpolation method or by projecting the local solution in each subdomain onto a global basis. In addition, we are working on extending the current approach to time-dependent PDEs. We will report on these issues in future publications.

214

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

Fig. 16. First three eigenfunctions of the covariance function of the solution for (27) with u (x, ξ ), d = 10 (in (a), (b), (c)) and with u (x, η), r = 3 (in (d), (e), (f)). (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

Acknowledgements This work was supported by the U.S. Department of Energy (DOE) Office of Advanced Scientific Computing. PNNL is operated by Battelle for the DOE under Contract DE-AC05-76RL01830. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19]

M. Loéve, Probability Theory, Springer-Verlag, 1977. R. Ghanem, P. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer-Verlag, 1991. D. Xiu, G.E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Sci. Comput. 24 (2002) 619–644. I. Babuška, P. Chatzipantelidis, On solving elliptic stochastic partial differential equations, Comput. Methods Appl. Mech. Eng. 191 (37–38) (2002) 4093–4122. I. Babuška, F. Nobile, R. Tempone, A stochastic collocation method for elliptic partial differential equations with random input data, SIAM Rev. 52 (2) (2010) 317–355. G. Lin, A. Tartakovsky, An efficient, high-order probabilistic collocation method on sparse grids for three-dimensional flow and solute transport in randomly heterogeneous porous media, Adv. Water Resour. 32 (5) (2009) 712–722. G. Lin, A.M. Tartakovsky, Numerical studies of three-dimensional stochastic Darcy’s equation and stochastic advection–diffusion–dispersion equation, J. Sci. Comput. 43 (1) (2010) 92–117. D. Venturi, D. Tartakovsky, A. Tartakovsky, G. Karniadakis, Exact PDF equations and closure approximations for advective–reactive transport, J. Comput. Phys. 243 (2013) 323–343. R. Tipireddy, R. Ghanem, Basis adaptation in homogeneous chaos spaces, J. Comput. Phys. 259 (2014) 304–317. P.G. Constantine, E. Dow, Q. Wang, Active subspace methods in theory and practice: applications to kriging surfaces, Methods Algorithms Sci. Comput. 36 (2014) A1500–A1524. A. Toselli, O.B. Widlund, Domain Decomposition Methods – Algorithms and Theory, Springer-Verlag, 2005. N. Agarwal, N. Aluru, A domain adaptive stochastic collocation approach for analysis of MEMS under uncertainties, J. Comput. Phys. 228 (20) (2009) 7662–7688. G. Lin, A. Tartakovsky, D. Tartakovsky, Uncertainty quantification via random domain decomposition and probabilistic collocation on sparse grids, J. Comput. Phys. 229 (19) (2010) 6995–7012. D. Xiu, D. Tartakovsky, A two-scale nonperturbative approach to uncertainty analysis of diffusion in random composites, Multiscale Model. Simul. 2 (4) (2004) 662–674. S. Acharjee, N. Zabaras, A concurrent model reduction approach on spatial and random domains for the solution of stochastic PDEs, Int. J. Numer. Methods Eng. 66 (12) (2006) 1934–1954. Y. Chen, J. Jakeman, C. Gittelson, D. Xiu, Local polynomial chaos expansion for linear differential equations with high dimensional random inputs, SIAM J. Sci. Comput. 37 (1) (2015) A79–A102. A. Doostan, R.G. Ghanem, J. Red-Horse, Stochastic model reduction for chaos representations, Comput. Methods Appl. Mech. Eng. 196 (2007) 3951–3966. R. Ghanem, The nonlinear Gaussian spectrum of log-normal stochastic processes and variables, J. Appl. Mech. 66 (4) (1999) 964. R.H. Cameron, W.T. Martin, The orthogonal development of non-linear functionals in series of Fourier–Hermite functionals, Ann. Math. 48 (2) (1947) 385.

R. Tipireddy et al. / Journal of Computational Physics 351 (2017) 203–215

[20] [21] [22] [23] [24]

215

A. Levy, J. Rubinstein, Hilbert-space Karhunen–Loéve transform with application to image analysis, J. Opt. Soc. Am. A 16 (1999) 28–35. M. Kirby, Minimal dynamical systems from PDEs using Sobolev eigenfunctions, Physica D 57 (1992) 466–475. B.W. Silverman, Smoothed functional principal components analysis by choice of norm, Ann. Stat. 24 (1996) 1–24. G. Berkooz, P. Holmes, J. Lumley, The proper orthogonal decomposition in the analysis of turbulent flows, Annu. Rev. Fluid Mech. 25 (1) (1993) 539–575. E.A. Christensen, M. Brøns, J.N. Sørensen, Evaluation of proper orthogonal decomposition–based decomposition techniques applied to parameterdependent nonturbulent flows, SIAM J. Sci. Comput. 21 (4) (1999) 1419–1434. [25] S. Smolyak, Quadrature and interpolation formulas for tensor products of certain classes of functions, Dokl. Akad. Nauk SSSR 4 (1963) 240–243. [26] F. Nobile, R. Tempone, C.G. Webster, A sparse grid stochastic collocation method for partial differential equations with random input data, SIAM J. Numer. Anal. 46 (2008) 2309–2345.