Stochastic Galerkin Methods for the Steady-State Navier-Stokes Equations This work is based upon work supported by the U.  S.  Department of Energy Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC0009301, and by the U.  S.  National Science Foundation under grants DMS1418754 and DMS1521563.

# Stochastic Galerkin Methods for the Steady-State Navier-Stokes Equations ††thanks: This work is based upon work supported by the U.  S.  Department of Energy Office of Advanced Scientific Computing Research, Applied Mathematics program under Award Number DE-SC0009301, and by the U.  S.  National Science Foundation under grants DMS1418754 and DMS1521563.

Bedřich Sousedík Department of Mathematics and Statistics, University of Maryland, Baltimore County, 1000 Hilltop Circle, Baltimore, MD 21250 (sousedik@umbc.edu).    Howard C. Elman Department of Computer Science and Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742 (elman@cs.umd.edu)
###### Abstract

We study the steady-state Navier-Stokes equations in the context of stochastic finite element discretizations. Specifically, we assume that the viscosity is a random field given in the form of a generalized polynomial chaos expansion. For the resulting stochastic problem, we formulate the model and linearization schemes using Picard and Newton iterations in the framework of the stochastic Galerkin method, and we explore properties of the resulting stochastic solutions. We also propose a preconditioner for solving the linear systems of equations arising at each step of the stochastic (Galerkin) nonlinear iteration and demonstrate its effectiveness for solving a set of benchmark problems.

## 1 Introduction

Models of mathematical physics are typically based on partial differential equations (PDEs) that use parameters as input data. In many situations, the values of parameters are not known precisely and are modeled as random fields, giving rise to stochastic partial differential equations. In this study we focus on models from fluid dynamics, in particular the stochastic Stokes and the Navier-Stokes equations. We consider the viscosity as a random field modeled as colored noise, and we use numerical methods based on spectral methods, specifically, the generalized polynomial chaos (gPC) framework [10, 14, 26, 27]. That is, the viscosity is given by a gPC expansion, and we seek gPC expansions of the velocity and pressure solutions.

There is a number of reasons to motivate our interest in Navier-Stokes equations with stochastic viscosity. For example, the exact value of viscosity may not be known, due to measurement error, the presence of contaminants with uncertain concentrations, or of multiple phases with uncertain ratios. Alternatively, the fluid properties might be influenced by an external field, with applications for example in magnetohydrodynamics. Specifically, we assume that the viscosity depends on a set of random variables . This means that the Reynolds number,

 Re(ξ)=ULν(ξ),

where is the characteristic velocity and is the characteristic length, is also stochastic. Consequently, the solution variables are random fields, and different realizations of the viscosity give rise to realizations of the velocities and pressures. As observed in [19], there are other possible formulations and interpretations of fluid flow with stochastic Reynolds number for example, where the velocity is fixed but the volume of fluid moving into a channel is uncertain so the uncertainty derives from the Dirichlet inflow boundary condition.

We consider models of steady-state stochastic motion of an incompressible fluid moving in a domain . Extension to three-dimensional models is straightforward. We formulate the stochastic Stokes and Navier-Stokes equations using the stochastic finite element method, assuming that the viscosity has a general probability distribution parametrized by a gPC expansion. We describe linearization schemes based on Picard and Newton iteration for the stochastic Galerkin method, and we explore properties of the solutions obtained, including a comparison of the stochastic Galerkin solutions with those obtained using other approaches, such as Monte Carlo and stochastic collocation methods [26]. Finally, we propose efficient hierarchical preconditioners for iterative solution of the linear systems solved at each step of the nonlinear iteration in the context of the stochastic Galerkin method. Our approach is related to recent work by Powell and Silvester [19]. However, besides using a general parametrization of the viscosity, our formulation of the stochastic Galerkin system allows straightforward application of state-of-the-art deterministic preconditioners by wrapping them in the hierarchical preconditioner developed in [22]. For alternative preconditioners see, e.g., [2, 11, 17, 18, 20, 23, 25]. Finally, we note that there exist related approaches based on stochastic perturbation methods [13], important developments also include reduced-order models such as [4, 24], and an overview of existing methods for stochastic computational fluid dynamics can be found in the monograph [14].

The paper is organized as follows. In Section LABEL:sec:NS, we recall the deterministic steady-state Navier-Stokes equations and their discrete form. In Section LABEL:sec:stochastic, we formulate the model with stochastic viscosity, derive linearization schemes for the stochastic Galerkin formulation of the model, and explore properties of the resulting solutions for a set of benchmark problems that model the flow over an obstacle. In Section LABEL:sec:preconditioner we introduce a preconditioner for the stochastic Galerkin linear systems solved at each step of the nonlinear iteration, and in Section LABEL:sec:conclusion we summarize our work.

## 2 Deterministic Navier-Stokes equations

We begin by defining the model and notation, following [6]. For the deterministic Navier-Stokes equations, we wish to find velocity and pressure such that

 −ν∇2→u+(→u⋅∇)→u+∇p =→f, \hb@xt@.01(2.1) ∇⋅→u =0, \hb@xt@.01(2.2)

in a spatial domain, satisfying boundary conditions

 →u =→g,on ΓDir, \hb@xt@.01(2.3) ν∇→u⋅→n−p→n =→0,on ΓNeu, \hb@xt@.01(2.4)

where, and assuming sufficient regularity of the data. Dropping the convective term from (LABEL:eq:NS-1) yields the Stokes problem

 −ν∇2→u+∇p =→f, \hb@xt@.01(2.5) ∇⋅→u =0. \hb@xt@.01(2.6)

The mixed variational formulation of (LABEL:eq:NS-1)–(LABEL:eq:NS-2) is to find such that

 ν∫D∇→u:∇→v+∫D(→u⋅∇→u)→v−∫Dp(∇⋅→v) =∫D→f⋅→v,∀→v∈VD, \hb@xt@.01(2.7) ∫Dq(∇⋅→u) =0,∀q∈QD, \hb@xt@.01(2.8)

where is a pair of spaces satisfying the inf-sup condition and is an extension of containing velocity vectors that satisfy the Dirichlet boundary conditions [3, 6, 12].

Let . Because the problem (LABEL:eq:NS-variational-1)–(LABEL:eq:NS-variational-2) is nonlinear, it is solved using a linearization scheme in the form of Newton or Picard iteration, derived as follows. Consider the solution of (LABEL:eq:NS-variational-1)–(LABEL:eq:NS-variational-2) to be given as and . Substituting into (LABEL:eq:NS-variational-1)–(LABEL:eq:NS-variational-2) and neglecting the quadratic term gives

 ν∫D∇δ→un:∇→v+c(δ→un;→un,→v)+c(→un;δ→un,→v)−∫Dδpn(∇⋅→v) =Rn(→v), \hb@xt@.01(2.9) ∫Dq(∇⋅δ→un) =rn(q), \hb@xt@.01(2.10)

where

 Rn(→v) =∫D→f⋅→v−ν∫D∇→un:∇→v−c(→un;→un,→v)+∫Dpn(∇⋅→v), \hb@xt@.01(2.11) rn(q) =−∫Dq(∇⋅→un). \hb@xt@.01(2.12)

Step  of the Newton iteration obtains from (LABEL:eq:Newton-1)–(LABEL:eq:Newton-2) and updates the solution as

 →un+1 =→un+δ→un, \hb@xt@.01(2.13) pn+1 =pn+δpn. \hb@xt@.01(2.14)

Step  of the Picard iteration omits the term in (LABEL:eq:Newton-1), giving

 ν∫D∇δ→un:∇→v+c(→un;δ→un,→v)−∫Dδpn(∇⋅→v) =Rn(→v), \hb@xt@.01(2.15) ∫Dq(∇⋅δ→un) =rn(q). \hb@xt@.01(2.16)

Consider the discretization of (LABEL:eq:NS-1)–(LABEL:eq:NS-2) by a div-stable mixed finite element method; for experiments discussed below, we used Taylor-Hood elements [6]. Let the bases for the velocity and pressure spaces be denoted  and , respectively. In matrix terminology, each nonlinear iteration entails solving a linear system

 [FnBTB0][δunδpn]=[Rnrn], \hb@xt@.01(2.17)

followed by an update of the solution

 un+1 =un+δun, \hb@xt@.01(2.18) pn+1 =pn+δpn. \hb@xt@.01(2.19)

For Newton’s method, is the Jacobian matrix, a sum of the vector-Laplacian matrix , the vector-convection matrix , and the Newton derivative matrix ,

 Fn=A+Nn+Wn, \hb@xt@.01(2.20)

where

 A=[aab], aab=ν∫D∇ϕb:∇ϕa, Nn=[nnab], nnab=∫D(un⋅∇ϕb)⋅ϕa, Wn=[wnab], wnab=∫D(ϕb⋅∇un)⋅ϕa.

For Picard iteration, the Newton derivative matrix  is dropped, and . The divergence matrix  is defined as

 B=[bcd],bcd=∫Dϕd(∇⋅φc). \hb@xt@.01(2.21)

The residuals at step  of both nonlinear iterations are computed as

 [Rnrn]=[fg]−[PnBTB0][unpn], \hb@xt@.01(2.22)

where and is a discrete version of the forcing function of (LABEL:eq:NS-1).111Throughout this study, we use the convention that the right-hand sides of discrete systems incorporate Dirichlet boundary data for velocities.

## 3 The Navier-Stokes equations with stochastic viscosity

Let represent a complete probability space, where is the sample space, is a-algebra on  and is a probability measure. We will assume that the randomness in the model is induced by a vector of independent, identically distributed (i.i.d.) random variables such that . Let  denote the -algebra generated by, and let denote the joint probability density measure for . The expected value of the product of random variables and that depend on determines a Hilbert space  with inner product

 ⟨u,v⟩=E[uv]=∫Γu(ξ)v(ξ)dμ(ξ), \hb@xt@.01(3.1)

where the symbol  denotes mathematical expectation.

### 3.1 The stochastic Galerkin formulation

The counterpart of the variational formulation (LABEL:eq:NS-variational-1)–(LABEL:eq:NS-variational-2) consists of performing a Galerkin projection on the space using mathematical expectation in the sense of (LABEL:eq:E). That is, we seek the velocity , a random field in , and the pressure , such that

 E[∫Dν∇→u:∇→v+∫D(→u⋅∇→u)→v−∫Dp(∇⋅→v)] =E[∫D→f⋅→v]∀→v∈VD⊗TΓ, \hb@xt@.01(3.2) E[∫Dq(∇⋅→u)] =0∀q∈QD⊗TΓ. \hb@xt@.01(3.3)

The stochastic counterpart of the Newton iteration (LABEL:eq:Newton-1)–(LABEL:eq:Newton-2) is

 E[∫Dν∇δ→un:∇→v+c(→un;δ→un,→v)+c(δ→un;→un,→v)−∫Dδpn(∇⋅→v)] =Rn, \hb@xt@.01(3.4) =rn, \hb@xt@.01(3.5)

where

 Rn(→v) =E[∫D→f⋅→v−∫Dν∇→un:∇→v−c(→un;→un,→v)+∫Dpn(∇⋅→v)], \hb@xt@.01(3.6) rn(q) =−E[∫Dq(∇⋅→un)]. \hb@xt@.01(3.7)

The analogue for Picard iteration omits from (LABEL:eq:Newton-1-s):

 E[∫Dν∇δ→un:∇→v+c(→un;δ→un,→v)−∫Dδpn(∇⋅→v)]=Rn. \hb@xt@.01(3.8)

In computations, we will use a finite-dimensional subspace spanned by a set of polynomials that are orthogonal with respect to the density function , that is . This is referred to as the gPC basis; see [10, 27] for details and discussion. For , we will use the space spanned by multivariate polynomials in  of total degree , which has dimension . We will also assume that the viscosity is given by a gPC expansion

 ν=Mν−1∑ℓ=0νℓ(x)ψℓ(ξ), \hb@xt@.01(3.9)

where is a set of given deterministic spatial functions.

### 3.2 Stochastic Galerkin finite element formulation

We discretize (LABEL:eq:Newton-1-s) (or (LABEL:eq:Picard-s)) and (LABEL:eq:Newton-2-s) using div-stable finite elements as in Section LABEL:sec:NS together with the gPC basis for . For simplicity, we assume that the right-hand side  and the Dirichlet boundary conditions (LABEL:eq:Dir-bc) are deterministic. This means in particular that, as in the deterministic case, the boundary conditions can be incorporated into right-hand side vectors (specified as  below). Thus, we seek a discrete approximation of the solution of the form

 →u(x,ξ) ≈M−1∑k=0Nu∑i=1uikϕi(x)ψk(ξ)=M−1∑k=0→uk(x)ψk(ξ), \hb@xt@.01(3.10) p(x,ξ) ≈M−1∑k=0Np∑j=1pjkφj(x)ψk(ξ)=M−1∑k=0pk(x)ψk(ξ), \hb@xt@.01(3.11)

The structure of the discrete operators depends on the ordering of the unknown coefficients , . We will group velocity-pressure pairs for each , the index of stochastic basis functions (and order equations in the same way), giving the ordered list of coefficients

 u1:Nu,0,p1:Np,0,u1:Nu,1,p1:Np,1,…,u1:Nu,M−1,p1:Np,M−1. \hb@xt@.01(3.12)

To describe the discrete structure, we first consider the stochastic version of the Stokes problem (LABEL:eq:S-1)–(LABEL:eq:S-2), where the convection term  is not present in (LABEL:eq:Newton-1-s) and (LABEL:eq:Picard-s). The discrete stochastic Stokes operator is built from the discrete components of the vector-Laplacian

 Aℓ=[aℓ,ab],aℓ,ab=(∫Dνℓ(x)∇ϕb:∇ϕa),ℓ=1,…,Mν−1, \hb@xt@.01(3.13)

which are incorporated into the block matrices

 \hb@xt@.01(3.14)

These operators will be coupled with matrices arising from terms in ,

 Hℓ=[hℓ,jk],hℓ,jk≡E[ψℓψjψk],ℓ=0,…,Mν−1,j,k=0,…,M−1. \hb@xt@.01(3.15)

Combining the expressions from (LABEL:stochastic-vector-Laplacian), (LABEL:stochastic-saddle-point) and (LABEL:eq:stochastic-matrix-forms) and using the ordering (LABEL:eq:coefficient-order) gives the discrete stochastic Stokes system

 (Mν−1∑ℓ=0Hℓ⊗Sℓ)v=y, \hb@xt@.01(3.16)

where corresponds to the matrix Kronecker product. The unknown vector corresponds to the ordered list of coefficients in (LABEL:eq:coefficient-order) and the right-hand side is ordered in an analogous way. Note that  is the identity matrix of order .

###### Remark 3.1

With this ordering, the coefficient matrix contains a set of  block  matrices of saddle-point structure along its block diagonal, given by

 S0+Mν−1∑ℓ=1hℓ,jjSℓ,j=0,…,M−1.

This enables the use of existing deterministic solvers for the individual diagonal blocks. An alternative ordering based on the blocking of all velocity coefficients followed by all pressure coefficients, considered in [19], produces a matrix of global saddle-point structure.

The matrices arising from the linearized stochastic Navier-Stokes equations augment the Stokes systems with stochastic variants of the vector-convection matrix and Newton derivative matrix appearing in (LABEL:Newton-Jacobian). In particular, at step  of the nonlinear iteration, let  be the th term of the velocity iterate (as in the expression on the right in (LABEL:eq:u-gPC) for ), and let

 Nnℓ=[nnℓ,ab], nnℓ,ab=∫D(→unℓ⋅∇ϕb)⋅ϕa, Wnℓ=[wnℓ,ab], wnℓ,ab=∫D(ϕb⋅∇→unℓ)⋅ϕa.

Then the analogues of (LABEL:stochastic-vector-Laplacian)–(LABEL:stochastic-saddle-point) are

 Fnℓ =Aℓ+Nnℓ+Wℓn, for the stochastic Newton method \hb@xt@.01(3.17) Fnℓ =Aℓ+Nnℓ, for stochastic Picard iteration, \hb@xt@.01(3.18)

so for Newton’s method

 Fn0=[Fn0BTB0],Fnℓ=[Fℓ000], \hb@xt@.01(3.19)

and as above, for Picard iteration the Newton derivative matrices  are dropped. Note that here, where . (In particular, if , we set for .) Step  of the stochastic nonlinear iteration entails solving a linear system and updating,

 ⎡⎣ˆM−1∑ℓ=0Hℓ⊗Fnℓ⎤⎦δvn=Rn,vn+1=vn+δvn, \hb@xt@.01(3.20)

where

 Rn=y−⎡⎣ˆM−1∑ℓ=0Hℓ⊗Pnℓ⎤⎦vn, \hb@xt@.01(3.21)

and are vectors of current velocity and pressure coefficients and updates, respectively, ordered as in (LABEL:eq:coefficient-order), is the similarly ordered right-hand side determined from the forcing function and Dirichlet boundary data, and

 Pn0=[A0+Nn0BTB0],Pnℓ=[Aℓ+Nnℓ000];

note that the -blocks here are as in (LABEL:eq:stoch-Picard).

### 3.3 Sampling methods

In experiments described below, we compare some results obtained using stochastic Galerkin methods to those obtained from Monte Carlo and stochastic collocation. We briefly describe these approaches here.

Both Monte Carlo and stochastic collocation methods are based on sampling. This entails the solution of a number of mutually independent deterministic problems at a set of sample points , which give realizations of the viscosity (LABEL:eq:viscosity-gPC). That is, a realization of viscosity  gives rise to deterministic functions  and on  that satisfy the standard deterministic Navier-Stokes equations, and to finite-element approximations .

In the Monte Carlo method, the sample points are generated randomly, following the distribution of the random variables, and moments of the solution are obtained from ensemble averaging. For stochastic collocation, the sample points consist of a set of predetermined collocation points. This approach derives from a methodology for performing quadrature or interpolation in multidimensional space using a small number of points, a so-called sparse grid [8, 16]. There are two ways to implement stochastic collocation to obtain the coefficients in (LABEL:eq:u-gPC)–(LABEL:eq:p-gPC), either by constructing a Lagrange interpolating polynomial, or, in the so-called pseudospectral approach, by performing a discrete projection into  [26]. We use the second approach because it facilitates a direct comparison with the stochastic Galerkin method. In particular, the coefficients are determined using a quadrature

 uik=Nq∑q=1→u(q)(xi)ψk(ξ(q))w(q),pik=Nq∑q=1p(q)(xi)ψk(ξ(q))w(q),

where  are collocation (quadrature) points, and  are quadrature weights. We refer, e.g., to [14] for an overview and discussion of integration rules.

### 3.4 Example: flow around an obstacle

In this section, we present results of numerical experiments for a model problem given by a flow around an obstacle in a channel of length and height . The viscosity (LABEL:eq:viscosity-gPC) was taken to be a truncated lognormal process with mean values or , which corresponds to mean Reynolds numbers or , respectively, and its representation was computed from an underlying Gaussian random process using the transformation described in [9]. That is, for , is the product of univariate Hermite polynomials, and denoting the coefficients of the Karhunen-Loève expansion of the Gaussian process by  and , , the coefficients in the expansion (LABEL:eq:viscosity-gPC) are computed as

 νℓ(x)=E[ψℓ(η)]exp[g0(x)+12N∑j=1(gj(x))2].

The covariance function of the Gaussian field, for points , , was chosen to be

 C(X1,X2)=σ2gexp(−|x2−x1|Lx−|y2−y1|Ly),

where  and  are the correlation lengths of the random variables , , in the and directions, respectively, and is the standard deviation of the Gaussian random field. The correlation lengths were set to be equal to of the width and height of the domain, i.e. and . The coefficient of variation of the lognormal field, defined as where is the standard deviation, was or . The stochastic dimension was . \textcolorblackThe lognormal formulation was chosen to enable exploration of a general random field in which the viscosity guaranteed to be positive. See [1] for an example of the use of this formulation for diffusion problems and [28] for its use in models of porous media.

We implemented the methods in Matlab using IFISS 3.3 [5]. The spatial discretization uses a stretched grid, discretized by Taylor-Hood finite elements; the domain and grid are shown in Figure LABEL:fig:mesh-obstacle. There are velocity and pressure degrees of freedom. the degree used for the polynomial expansion of the solution was , and the degree used for the expansion of the lognormal process was , which ensures a complete representation of the process in the discrete problem [15]. With these settings, and , and is of order in (LABEL:eq:linearized-system).

Consider first the case of and . Figure LABEL:fig:sRe100_CoV10_0 shows the mean horizontal and vertical components of the velocity and the mean pressure (top), and the variances of the same quantities (bottom). It can be seen that there is symmetry in all the quantities, the mean values are essentially the same as we would expect in the deterministic case, and the variance of the horizontal velocity component is concentrated in two “eddies” and is larger than the variance of the vertical velocity component. Figure LABEL:fig:sRe100_CoV10_gPC illustrates values of several coefficients of expansion (LABEL:eq:u-gPC) of the horizontal velocity. All the coefficients are symmetric, and as the index increases they become more oscillatory and their values decay. We found the same trends for the coefficients of the vertical velocity component and of the pressure. Our observations are qualitatively consistent with numerical experiments of Powell and Silvester [19].

We also tested (the same) with increased . We found that the mean values are essentially the same as in the previous case; Figure LABEL:fig:sRe100_CoV30_var shows the variances, which display the same qualitative behavior but have values that are approximately times larger than for the case .

A different perspective on the results is given in Figure LABEL:fig:sRe100_QoI, which shows estimates of the probability density function (pdf) for the horizontal velocity at two points in the domain, and . These are locations at which large variances of the solution were seen, see Figures LABEL:fig:sRe100_CoV10_0 and LABEL:fig:sRe100_CoV30_var. The results were obtained using Matlab’s ksdensity function. It can be seen that with the larger value of , the support of the velocity pdf is wider, and except for the peak values, for fixed the shapes of the pdfs at the two points are similar, indicating a possible symmetry of the stochastic solution. For this benchmark, we also obtained analogous data using the Monte Carlo and collocation sampling methods; it can be seen from the figure that these methods produced similar results. 222The results for Monte Carlo were obtained using samples, and those for collocation were found using a Smolyak sparse grid with Gauss-Hermite quadrature and grid level .

Next, we consider a larger value of the mean Reynolds number, . Figure LABEL:fig:sRe300_CoV10_0 shows the means and variances for the velocities and pressure for . It is evident that increased results in increased values of the mean quantities, but they are again similar to what would be expected in the deterministic case. The variances exhibit wider eddies than for , and in this case there is only one region of the largest variance in the horizontal velocity, located just to the right of the obstacle; this is also a region with increased variance of the pressure.

In similar tests for the larger value , we found that the mean values are essentially the same as for , and Figure LABEL:fig:sRe300_CoV30_var shows the variances of velocities and pressures. From the figure it can be seen that the variances are qualitatively the same but approximately times larger than for , results similar to those found for .

As above, we also examined estimated probability density functions for the velocities and pressures at a specified point in the domain, in this case, at the point , taken again from the region in which the solution has large variance. Figure LABEL:fig:sRe300_QoI shows these pdf estimates, from which it can be seen that the three methods for handling uncertainty are in close agreement for each of the two values of .

\textcolor

blackFinally, we show in Figure LABEL:fig:inflow_QoI the results of one additional experiment with and , where estimated pdfs of the first velocity component  were computed at three points near the inflow boundary, , and . These plots show some effects of spatial accuracy. \colorblackThe results for all methods were identical, and so only one representative pdf is shown. The images on the top and bottom left show results for a uniform mesh of width and for two refined meshes in which the horizontal mesh width to the left of the obstacle is reduced to and . The image in the bottom right provides a more detailed view of the fine-grid results; in this image, the width of the horizontal window is the same () for the three subplots but the vertical heights are different. Several trends are evident:

• The pdfs at points nearer the boundary are much narrower.

• When the spatial accuracy is low, \colorblack the methods miss some features of the pdf. \colorblackIn particular the methods produce a pdf that captures its “spiky” nature near the inflow, but the location of the spike is not correct.

We believe these effects stem from the fact that the inflow boundary is deterministic (with at ), and the effects of variability in the viscosity are felt less strongly at points near the inflow boundary. At points further from the deterministic inflow boundary, these effects and differences become less dramatic.

### 3.5 Nonlinear solvers

We briefly comment on the nonlinear solution algorithm used to generate the results of the previous section. The nonlinear solver was implemented by modifying the analogue for deterministic systems in IFISS. It uses a hybrid strategy in which an initial approximation is obtained from solution of the stochastic Stokes problem (LABEL:eq:Stokes-s), after which several steps of Picard iteration (equation (LABEL:eq:linearized-system) with specified using (LABEL:eq:linearized-operator) and (LABEL:eq:stoch-Picard)) are used to improve the solution, followed by Newton iteration ( from (LABEL:eq:stoch-Newton)). A convergent iteration stopped when the Euclidian norm of the algebraic residual (LABEL:eq:R_gPC) satisfied where and is as in (LABEL:eq:Stokes-s).

In the experiments described in Section LABEL:sec:obstacle, we used values of the Reynolds number, and , and for each of these, two values of the coefficient of variation, and . We list here the numbers of steps leading to convergence of the nonlinear algorithms that were used to generate the solutions discussed above. Direct solvers were used for the linear systems; we discuss a preconditioned iterative algorithm in Section LABEL:sec:preconditioner below. , : Picard steps Newton step(s) , : 6 \colorblack2 , : 20 1 , : 20 2 Thus, a larger (larger standard deviation of the random field determining uncertainty in the process) leads to somewhat larger computational costs. For , the nonlinear iteration was not robust with initial Picard steps (for the stochastic Galerkin method as well as the sampling methods); steps was sufficient.

We also explored an inexact variant of these methods, in which the coefficient matrix of (LABEL:eq:linearized-system) for the Picard iteration was replaced by the block diagonal matrix obtained from the mean coefficient. For , with the same number of (now inexact) Picard steps as above ( for and for ), this led to just one extra (exact) Newton step for and no additional steps for . On the other hand, for , this inexact method failed to converge.

## 4 Preconditioner for the linearized systems

The solution of the linear systems required during the course of the nonlinear iteration is a computationally intensive task, and use of direct solvers may be prohibitive for large problems. In this section, we present a preconditioning strategy for use with Krylov subspace methods to solve these systems, and we compare its performance with that of several other techniques. The new method is a variant of the hierarchical Gauss-Seidel preconditioner developed in [22].

### 4.1 Structure of the matrices and the preconditioner

We first recall the structure of the matrices of (LABEL:eq:stochastic-matrix-forms). More comprehensive overviews of these matrices can be found in [7, 15]. The matrix structure can be understood through knowledge of the coefficient matrix where . The block sparsity structure depends on the type of coefficient expansion in (LABEL:eq:viscosity-gPC). If only linear terms are included, that is , , then the coefficients yield a Galerkin matrix with a block sparse structure. In the more general case, and the stochastic Galerkin matrix becomes fully block dense. In either case, for fixed and a set of degree polynomial expansions, with , the corresponding coefficient matrix  has a hierarchical structure

 cP=[cP−1bTPbPdP],P=1,…,P.

Now, let denote the global stochastic Galerkin matrix corresponding to either a Stokes problem (LABEL:eq:Stokes-s) or a linearized system (LABEL:eq:linearized-system); we will focus on the latter system in the discussion below. The matrix  also has a hierarchical structure

 AP=[AP−1BPCPDP],P=1,…,P, \hb@xt@.01(4.1)

where  is the matrix of the mean, derived from  in (LABEL:eq:viscosity-gPC). This hierarchical structure is shown in the left side of Figure LABEL:fig:GS.

We will write vectors with respect to this hierarchy as

 x(0:P)=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣x(0)x(1)⋮x(P)⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

where includes all indices corresponding to polynomial degree , blocked by spatial ordering determined by (LABEL:eq:coefficient-order). With this notation, the global stochastic Galerkin linear system has the form

 APx(0:P)=f(0:P). \hb@xt@.01(4.2)

To formulate the preconditioner for (LABEL:eq:SG), we let represent an approximation of  and represent an approximation of . In particular, let

 ˜DP=⎡⎢ ⎢ ⎢ ⎢⎣˜A0⋱˜A0⎤⎥ ⎥ ⎥ ⎥⎦, \hb@xt@.01(4.3)

where the number of diagonal blocks is given by . We will need the action of the inverse of , or an approximation to it, which can be obtained using an LU-factorization of , or using some preconditioner for , or using a Krylov subspace solver. A preconditioner  for (LABEL:eq:SG) is then defined as follows:

###### Algorithm 4.1

[Approximate hierarchical Gauss-Seidel preconditioner (ahGS)] Solve (or solve approximately)

 ˜A0v(0)=w(0), \hb@xt@.01(4.4)

and, for , solve (or solve approximately)

 ˜DPv(P)=(w(P)−CPv(0:P−1)). \hb@xt@.01(4.5)

The cost of preconditioning can be reduced further by truncating the matrix-vector (MATVEC) operations used for the multiplications by the submatrices  in (LABEL:eq:ahGS-2). The idea is as follows. The system (LABEL:eq:SG) can be written as

 M−1∑j=0Mν−1∑ℓ=0hℓ,jkFℓxj=fk,k=0,…,M−1, \hb@xt@.01(4.6)

and the MATVEC with is given by

 vj=M−1∑k=0Mν−1∑ℓ=0hℓ,jkFℓuk, \hb@xt@.01(4.7)

where the indices , correspond to nonzero blocks in. The truncated MATVEC is an inexact evaluation of (LABEL:eq:MAT-VEC) proposed in [22, Algorithm 1], in which the summation over is replaced by summation over a subset . Figure LABEL:fig:GS shows the hierarchical structure of the matrix and of the ahGS preconditioning operator. Both images in the figure correspond to the choice , so that the hierarchical preconditioning operation (LABEL:eq:ahGS-1)–(LABEL:eq:ahGS-2) requires four steps. Because , the matrix block size is . The block-lower-triangular component of the image on the right in Figure LABEL:fig:GS shows the hierarchical structure of the ahGS preconditioning operator with truncation. For the matrix in the left panel, , but the index set includes terms with indices at most in the accumulation of sums used for. These two heuristics, approximation by (LABEL:eq:D-approx) and the truncation of MATVECs, significantly improve the sparsity structure of the preconditioner, on both the block diagonal (through the first technique) and the block lower triangle (through the second). In the next section, we will also consider truncated MATVECs with smaller maximal indices.

### 4.2 Numerical experiments

In this section, we describe the results of experiments in which the ahGS preconditioner is used with GMRES to solve the systems arising from both Picard and Newton iteration.333The preconditioning operator may be nonlinear, for example, if the block solves in (LABEL:eq:ahGS-1)–(LABEL:eq:ahGS-2) are performed approximately using Krylov subspace methods, so that care must be made in the choice of the Krylov subspace iterative method. For this reason, we used a flexible variant of the GMRES method [21]. \textcolorblackUnless otherwise specified, in the experiments discussed below, the block-diagonal matrices of (LABEL:eq:D-approx), (LABEL:eq:ahGS-1), (LABEL:eq:ahGS-2) are taken to be , which corresponds to the first summand in (LABEL:eq:linearized-system) and represents an approximation to the diagonal blocks of ; see also Remark LABEL:rem:order. We also compare the performance of ahGS preconditioning with several alternatives:

• Mean-based preconditioner (MB) [17, 18], where the preconditioning operator is the block-diagonal matrix derived from the mean of .

• The Kronecker-product preconditioner [25] (denoted K), given by , where is chosen to to minimize a measure of the difference .

• Block Gauss-Seidel preconditioner (bGS), in which the preconditioning operation entails applying determined using the block lower triangle of , that is with the approximation of the diagonal blocks by but without MATVEC truncation. This is an expensive operator but enables an understanding of the impact of various efforts to make the preconditioning operator more sparse.

• \textcolor

blackahGS(PCD), a modification of ahGS in which the block-diagonal matrices of (LABEL:eq:ahGS-1)–(LABEL:eq:ahGS-2) are replaced by the pressure convection-diffusion (PCD) approximation [6] to the mean matrix .

• \textcolor

blackahGS(PCD-it), in which the block-diagonal solve of (LABEL:eq:ahGS-2) is replaced by an approximate solve determined by twenty steps of a PCD-preconditioned GMRES iteration. For this strategy, in (LABEL:eq:ahGS-1) and the blocks of in (LABEL:eq:ahGS-2) are taken to be the complete sum of (LABEL:eq:linearized-system), but the approximate solutions to these systems are obtained using a preconditioned inner iteration.

Four of the strategies, the ahGS, mean-based, Kronecker-product and bGS preconditioners, require the solution of a set of block-diagonal systems with the structure of a linearized Navier-Stokes operator of the form given in (LABEL:eq:Newton). We used direct methods for these computations. All results presented are for ; performance for were essentially the same. We believe this is because of the exact solves performed for the mean operators.

The results for the first step of Picard iteration are in Tables LABEL:tab:Re_100-Picard-NLABEL:tab:Re_100-Picard-trunc. All tests started with a zero initial iterate and stoppped when the residual for the ’th iterate satisfied in the Euclidian norm. With other parameters fixed and no truncation of the MATVEC, Table LABEL:tab:Re_100-Picard-N shows the dependence of GMRES iterations on the stochastic dimension , Table LABEL:tab:Re_100-Picard-P shows the dependence on the degree of polynomial expansion , and Table LABEL:tab:Re_100-Picard-CoV shows the dependence on the coefficient of variation . It can be seen that the numbers of iterations with the ahGS preconditioner are essentially the same as with the block Gauss-Seidel (bGS) preconditioner, and they are much smaller compared to the mean-based (MB) and the Kronecker product preconditioners. When the exact solves with the mean matrix are replaced by the mean-based modified pressure-convection-diffusion (PCD) preconditioner for the diagonal block solves, the iterations grow rapidly. \textcolorblackOn the other hand, if the PCD preconditioner is used as part of an inner iteration as in ahGS(PCD-it), then good performance is recovered. This indicates that a good preconditioner for the mean matrix is an essential component of the global preconditioner for the stochastic Galerkin matrix.

Table LABEL:tab:Re_100-Picard-trunc shows the iteration counts when the MATVEC operation is truncated in the action of the preconditioner. Truncation decreases the cost per iteration of the computation, and it can be also seen that performance can actually be improved. For example, with  , the number of iterations is the smallest. Moreover there are only nonzeros in the lower triangular part of the sum of coefficient matrices  (each of which has order ) used in the MATVEC with , compared to nonzeros when no truncation is used; there are nonzeros in the set of full matrices .

Results for the first step of Newton iteration, which comes after six steps of Picard iteration, are summarized in Tables LABEL:tab:Re_100-Newton-NLABEL:tab:Re_100-Newton-trunc. As above, the first three tables show the dependence on stochastic dimension (Table LABEL:tab:Re_100-Newton-N), polynomial degree (Table LABEL:tab:Re_100-Newton-P), and coefficient of variation (Table LABEL:tab:Re_100-Newton-CoV). It can be seen that all iteration counts are higher compared to corresponding results for the Picard iteration, but the other trends are very similar. In particular, the performances of the ahGS and bGS preconditioners are comparable except for the case when , and (the last rows in Tables LABEL:tab:Re_100-Newton-N and LABEL:tab:Re_100-Newton-P). Nevertheless, checking results in Table LABEL:tab:Re_100-Newton-trunc, which shows the effect of truncation, it can be seen that with the truncation of the MATVEC the iteration counts of the ahGS and bGS preconditioners can be further improved. Indeed, running one more experiment for the aforementioned case when , and , it turns out that with the number of iterations with the ahGS and bGS preconditioners are and , respectively and with they are and , respectively. That is, the truncation leads to fewer iterations also in this case, and the performance of ahGS and bGS preconditioners is again comparable.

Finally, we briefly discuss computational costs. For any preconditioner, each GMRES step entails a matrix-vector product by the coefficient matrix. For viscosity given by a general probability distribution, this will typically involve a block-dense matrix, and, ignoring any overhead associated with increasing the number of GMRES steps, this will be the dominant cost per step. The mean-based preconditioner requires the action of the inverse of the block-diagonal matrix . This has relatively small amount of overhead once the factors of are computed. The ahGS preconditioner without truncation effectively entails a matrix-vector product by the block lower-triangular part of the coefficient matrix, so its overhead is bounded by of the cost of a multiplication by the coefficient matrix. This is an overestimate because it ignores the approximation of the block diagonal and the effect of truncation. For example, consider the case with stochastic dimension and polynomial expansions of degree for the solution and for the viscosity; this gives and . In Tables LABEL:tab:Re_100-Picard-trunc and LABEL:tab:Re_100-Newton-trunc,  indicates how many matrices are used in the MATVEC operations and is the number of nonzeros in the sum of the lower triangular parts of the coefficient matrices . With complete truncation, , and ahGS reduces to the mean-based preconditioner. With no truncation, , and because the number of nonzeros in is , the overhead of ahGS is , less than of the cost of multiplication by the coefficient matrix. If truncation is used, in particular when the iteration count is the lowest (), the overhead is only . Note that with increasing stochastic dimension and degree of polynomial expansion, the savings will be higher because the ratio of the sizes of the blocks decreases as increases, see (LABEL:eq:hSG). Last, the mean-based preconditioner is embarrassingly parallel; the ahGS preconditioner requires sequential steps, although each of these steps is also highly parallelizable. The Kronecker preconditioner is more difficult to assess because it does not have block-diagonal structure, and we do not discuss it here.

## 5 Conclusion

We studied the Navier-Stokes equations with stochastic viscosity given in terms of polynomial chaos expansion. We formulated the stochastic Galerkin method and proposed its numerical solution using a stochastic versions of Picard and Newton iteration, and we also compared its performance in terms of accuracy with that of stochastic collocation and Monte Carlo method. Finally, we presented a methodology of Gauss-Seidel hierarchical preconditioning with approximation using the mean-based diagonal block solves and a truncation of the MATVEC operations. The advantage of this approach is that neither the matrix nor the preconditioner need to be formed explicitly, and the ingredients include only the matrices from the polynomial chaos expansion and a good preconditioner for the mean-value deterministic problem, it allows an obvious parallel implementation, and it can be written as a “wrapper” around existing deterministic code.

## References

• [1] I. Babuška, F. Nobile, and R. Tempone. A stochastic collocation method for elliptic partial differential equations with random input data. SIAM Journal on Numerical Analysis, 45:1005–1034, 2007.
• [2] M. Brezina, A. Doostan, T. Manteuffel, S. McCormick, and J. Ruge. Smoothed aggregation algebraic multigrid for stochastic PDE problems with layered materials. Numerical Linear Algebra with Applications, 21(2):239–255, 2014.
• [3] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods. Springer-Verlag, New York – Berlin – Heidelberg, 1991.
• [4] H. C. Elman and Q. Liao. Reduced basis collocation methods for partial differential equations with random coefficients. SIAM/ASA Journal on Uncertainty Quantification, 1:192–217, 2013.
• [5] H. C. Elman, A. Ramage, and D. J. Silvester. IFISS: A computational laboratory for investigating incompressible flow problems. SIAM Review, 56:261–273, 2014.
• [6] H. C. Elman, D. J. Silvester, and A. J. Wathen. Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Numerical Mathematics and Scientific Computation. Oxford University Press, New York, second edition, 2014.
• [7] O. G. Ernst and E. Ullmann. Stochastic Galerkin matrices. SIAM Journal on Matrix Analysis and Applications, 31(4):1848–1872, 2010.
• [8] T. Gerstner and M. Griebel. Numerical integration using sparse grids. Numerical Algorithms, 18:209–232, 1998.
• [9] R. Ghanem. The nonlinear Gaussian spectrum of log-normal stochastic processes and variables. Journal of Applied Mechanics, 66(4):964–973, 1999.
• [10] R. G. Ghanem and Pol D. Spanos. Stochastic Finite Elements: A Spectral Approach. Springer-Verlag New York, Inc., New York, NY, USA, 1991. (Revised edition by Dover Publications, 2003).
• [11] L. Giraldi, A. Nouy, and G. Legrain. Low-rank approximate inverse for preconditioning tensor-structured linear systems. SIAM Journal on Scientific Computing, 36(4):A1850–A1870, 2014.
• [12] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations. Springer-Verlag, Berlin, 1986.
• [13] M. Kamiński. The Stochastic Perturbation Method for Computational Mechanics. John Wiley & Sons, 2013.
• [14] O. Le Maître and O. M. Knio. Spectral Methods for Uncertainty Quantification: With Applications to Computational Fluid Dynamics. Scientific Computation. Springer, 2010.
• [15] H. G. Matthies and A. Keese. Galerkin methods for linear and nonlinear elliptic stochastic partial differential equations. Computer Methods in Applied Mechanics and Engineering, 194(12â16):1295–1331, 2005.
• [16] E. Novak and K. Ritter. High dimensional integration of smooth functions over cubes. Numerische Mathematik, 75(1):79–97, 1996.
• [17] M. F. Pellissetti and R. G. Ghanem. Iterative solution of systems of linear equations arising in the context of stochastic finite elements. Advances in Engineering Software, 31(8-9):607–616, 2000.
• [18] C. E. Powell and H C. Elman. Block-diagonal preconditioning for spectral stochastic finite-element systems. IMA Journal of Numerical Analysis, 29(2):350–375, 2009.
• [19] C. E. Powell and D. J. Silvester. Preconditioning steady-state Navier–Stokes equations with random data. SIAM Journal on Scientific Computing, 34(5):A2482–A2506, 2012.
• [20] E. Rosseel and S. Vandewalle. Iterative solvers for the stochastic finite element method. SIAM Journal on Scientific Computing, 32(1):372–397, 2010.
• [21] Y. Saad. A flexible inner-outer preconditioned GMRES algorithm. SIAM Journal on Scientific Computing, 14(2):461–469, 1993.
• [22] B. Sousedík and R. G. Ghanem. Truncated hierarchical preconditioning for the stochastic Galerkin FEM. International Journal for Uncertainty Quantification, 4(4):333–348, 2014.
• [23] B. Sousedík, R. G. Ghanem, and E. T. Phipps. Hierarchical Schur complement preconditioner for the stochastic Galerkin finite element methods. Numerical Linear Algebra with Applications, 21(1):136–151, 2014.
• [24] L. Tamellini, O. Le Maître, and A. Nouy. Model reduction based on Proper Generalized Decomposition for the stochastic steady incompressible Navier–Stokes equations. SIAM Journal on Scientific Computing, 36(3):A1089–A1117, 2014.
• [25] E. Ullmann. A Kronecker product preconditioner for stochastic Galerkin finite element discretizations. SIAM Journal on Scientific Computing, 32(2):923–946, 2010.