Stochastic Galerkin Methods for the SteadyState NavierStokes 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 DESC0009301, and by the U. S. National Science Foundation under grants DMS1418754 and DMS1521563.
Abstract
We study the steadystate NavierStokes 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 NavierStokes 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 NavierStokes 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,
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 steadystate stochastic motion of an incompressible fluid moving in a domain . Extension to threedimensional models is straightforward. We formulate the stochastic Stokes and NavierStokes 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 stateoftheart 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 reducedorder 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 steadystate NavierStokes 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 NavierStokes equations
We begin by defining the model and notation, following [6]. For the deterministic NavierStokes equations, we wish to find velocity and pressure such that
\hb@xt@.01(2.1)  
\hb@xt@.01(2.2) 
in a spatial domain, satisfying boundary conditions
\hb@xt@.01(2.3)  
\hb@xt@.01(2.4) 
where, and assuming sufficient regularity of the data. Dropping the convective term from (LABEL:eq:NS1) yields the Stokes problem
\hb@xt@.01(2.5)  
\hb@xt@.01(2.6) 
The mixed variational formulation of (LABEL:eq:NS1)–(LABEL:eq:NS2) is to find such that
\hb@xt@.01(2.7)  
\hb@xt@.01(2.8) 
where is a pair of spaces satisfying the infsup condition and is an extension of containing velocity vectors that satisfy the Dirichlet boundary conditions [3, 6, 12].
Let . Because the problem (LABEL:eq:NSvariational1)–(LABEL:eq:NSvariational2) 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:NSvariational1)–(LABEL:eq:NSvariational2) to be given as and . Substituting into (LABEL:eq:NSvariational1)–(LABEL:eq:NSvariational2) and neglecting the quadratic term gives
\hb@xt@.01(2.9)  
\hb@xt@.01(2.10) 
where
\hb@xt@.01(2.11)  
\hb@xt@.01(2.12) 
Step of the Newton iteration obtains from (LABEL:eq:Newton1)–(LABEL:eq:Newton2) and updates the solution as
\hb@xt@.01(2.13)  
\hb@xt@.01(2.14) 
Step of the Picard iteration omits the term in (LABEL:eq:Newton1), giving
\hb@xt@.01(2.15)  
\hb@xt@.01(2.16) 
Consider the discretization of (LABEL:eq:NS1)–(LABEL:eq:NS2) by a divstable mixed finite element method; for experiments discussed below, we used TaylorHood 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
\hb@xt@.01(2.17) 
followed by an update of the solution
\hb@xt@.01(2.18)  
\hb@xt@.01(2.19) 
For Newton’s method, is the Jacobian matrix, a sum of the vectorLaplacian matrix , the vectorconvection matrix , and the Newton derivative matrix ,
\hb@xt@.01(2.20) 
where
For Picard iteration, the Newton derivative matrix is dropped, and . The divergence matrix is defined as
\hb@xt@.01(2.21) 
The residuals at step of both nonlinear iterations are computed as
\hb@xt@.01(2.22) 
where and is a discrete version of the forcing function of (LABEL:eq:NS1).^{1}^{1}1Throughout this study, we use the convention that the righthand sides of discrete systems incorporate Dirichlet boundary data for velocities.
3 The NavierStokes equations with stochastic viscosity
Let represent a complete probability space, where is the sample space, is aalgebra 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
\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:NSvariational1)–(LABEL:eq:NSvariational2) 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
\hb@xt@.01(3.2)  
\hb@xt@.01(3.3) 
The stochastic counterpart of the Newton iteration (LABEL:eq:Newton1)–(LABEL:eq:Newton2) is
\hb@xt@.01(3.4)  
\hb@xt@.01(3.5) 
where
\hb@xt@.01(3.6)  
\hb@xt@.01(3.7) 
The analogue for Picard iteration omits from (LABEL:eq:Newton1s):
\hb@xt@.01(3.8) 
In computations, we will use a finitedimensional 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
\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:Newton1s) (or (LABEL:eq:Picards)) and (LABEL:eq:Newton2s) using divstable finite elements as in Section LABEL:sec:NS together with the gPC basis for . For simplicity, we assume that the righthand side and the Dirichlet boundary conditions (LABEL:eq:Dirbc) are deterministic. This means in particular that, as in the deterministic case, the boundary conditions can be incorporated into righthand side vectors (specified as below). Thus, we seek a discrete approximation of the solution of the form
\hb@xt@.01(3.10)  
\hb@xt@.01(3.11) 
The structure of the discrete operators depends on the ordering of the unknown coefficients , . We will group velocitypressure pairs for each , the index of stochastic basis functions (and order equations in the same way), giving the ordered list of coefficients
\hb@xt@.01(3.12) 
To describe the discrete structure, we first consider the stochastic version of the Stokes problem (LABEL:eq:S1)–(LABEL:eq:S2), where the convection term is not present in (LABEL:eq:Newton1s) and (LABEL:eq:Picards). The discrete stochastic Stokes operator is built from the discrete components of the vectorLaplacian
\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 ,
\hb@xt@.01(3.15) 
Combining the expressions from (LABEL:stochasticvectorLaplacian), (LABEL:stochasticsaddlepoint) and (LABEL:eq:stochasticmatrixforms) and using the ordering (LABEL:eq:coefficientorder) gives the discrete stochastic Stokes system
\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:coefficientorder) and the righthand 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 saddlepoint structure along its block diagonal, given by
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 saddlepoint structure.
The matrices arising from the linearized stochastic NavierStokes equations augment the Stokes systems with stochastic variants of the vectorconvection matrix and Newton derivative matrix appearing in (LABEL:NewtonJacobian). 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:ugPC) for ), and let
Then the analogues of (LABEL:stochasticvectorLaplacian)–(LABEL:stochasticsaddlepoint) are
for the stochastic Newton method  \hb@xt@.01(3.17)  
for stochastic Picard iteration,  \hb@xt@.01(3.18) 
so for Newton’s method
\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,
\hb@xt@.01(3.20) 
where
\hb@xt@.01(3.21) 
and are vectors of current velocity and pressure coefficients and updates, respectively, ordered as in (LABEL:eq:coefficientorder), is the similarly ordered righthand side determined from the forcing function and Dirichlet boundary data, and
note that the blocks here are as in (LABEL:eq:stochPicard).
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:viscositygPC). That is, a realization of viscosity gives rise to deterministic functions and on that satisfy the standard deterministic NavierStokes equations, and to finiteelement 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 socalled sparse grid [8, 16]. There are two ways to implement stochastic collocation to obtain the coefficients in (LABEL:eq:ugPC)–(LABEL:eq:pgPC), either by constructing a Lagrange interpolating polynomial, or, in the socalled 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
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:viscositygPC) 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 KarhunenLoève expansion of the Gaussian process by and , , the coefficients in the expansion (LABEL:eq:viscositygPC) are computed as
The covariance function of the Gaussian field, for points , , was chosen to be
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 TaylorHood finite elements; the domain and grid are shown in Figure LABEL:fig:meshobstacle. 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:linearizedsystem).
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:ugPC) 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. ^{2}^{2}2The results for Monte Carlo were obtained using samples, and those for collocation were found using a Smolyak sparse grid with GaussHermite 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 .
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 finegrid 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:Stokess), after which several steps of Picard iteration (equation (LABEL:eq:linearizedsystem) with specified using (LABEL:eq:linearizedoperator) and (LABEL:eq:stochPicard)) are used to improve the solution, followed by Newton iteration ( from (LABEL:eq:stochNewton)). A convergent iteration stopped when the Euclidian norm of the algebraic residual (LABEL:eq:R_gPC) satisfied where and is as in (LABEL:eq:Stokess).
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:linearizedsystem) 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 GaussSeidel preconditioner developed in [22].
4.1 Structure of the matrices and the preconditioner
We first recall the structure of the matrices of (LABEL:eq:stochasticmatrixforms). 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:viscositygPC). 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
Now, let denote the global stochastic Galerkin matrix corresponding to either a Stokes problem (LABEL:eq:Stokess) or a linearized system (LABEL:eq:linearizedsystem); we will focus on the latter system in the discussion below. The matrix also has a hierarchical structure
\hb@xt@.01(4.1) 
where is the matrix of the mean, derived from in (LABEL:eq:viscositygPC). This hierarchical structure is shown in the left side of Figure LABEL:fig:GS.
We will write vectors with respect to this hierarchy as
where includes all indices corresponding to polynomial degree , blocked by spatial ordering determined by (LABEL:eq:coefficientorder). With this notation, the global stochastic Galerkin linear system has the form
\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
\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 LUfactorization 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 GaussSeidel preconditioner (ahGS)] Solve (or solve approximately)
\hb@xt@.01(4.4) 
and, for , solve (or solve approximately)
\hb@xt@.01(4.5) 
The cost of preconditioning can be reduced further by truncating the matrixvector (MATVEC) operations used for the multiplications by the submatrices in (LABEL:eq:ahGS2). The idea is as follows. The system (LABEL:eq:SG) can be written as
\hb@xt@.01(4.6) 
and the MATVEC with is given by
\hb@xt@.01(4.7) 
where the indices , correspond to nonzero blocks in. The truncated MATVEC is an inexact evaluation of (LABEL:eq:MATVEC) 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:ahGS1)–(LABEL:eq:ahGS2) requires four steps. Because , the matrix block size is . The blocklowertriangular 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:Dapprox) 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.^{3}^{3}3The preconditioning operator may be nonlinear, for example, if the block solves in (LABEL:eq:ahGS1)–(LABEL:eq:ahGS2) 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 blockdiagonal matrices of (LABEL:eq:Dapprox), (LABEL:eq:ahGS1), (LABEL:eq:ahGS2) are taken to be , which corresponds to the first summand in (LABEL:eq:linearizedsystem) 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:

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

Block GaussSeidel 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 blockdiagonal matrices of (LABEL:eq:ahGS1)–(LABEL:eq:ahGS2) are replaced by the pressure convectiondiffusion (PCD) approximation [6] to the mean matrix .

\textcolor
blackahGS(PCDit), in which the blockdiagonal solve of (LABEL:eq:ahGS2) is replaced by an approximate solve determined by twenty steps of a PCDpreconditioned GMRES iteration. For this strategy, in (LABEL:eq:ahGS1) and the blocks of in (LABEL:eq:ahGS2) are taken to be the complete sum of (LABEL:eq:linearizedsystem), but the approximate solutions to these systems are obtained using a preconditioned inner iteration.
Four of the strategies, the ahGS, meanbased, Kroneckerproduct and bGS preconditioners, require the solution of a set of blockdiagonal systems with the structure of a linearized NavierStokes 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_100PicardN–LABEL:tab:Re_100Picardtrunc. 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_100PicardN shows the dependence of GMRES iterations on the stochastic dimension , Table LABEL:tab:Re_100PicardP shows the dependence on the degree of polynomial expansion , and Table LABEL:tab:Re_100PicardCoV 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 GaussSeidel (bGS) preconditioner, and they are much smaller compared to the meanbased (MB) and the Kronecker product preconditioners. When the exact solves with the mean matrix are replaced by the meanbased modified pressureconvectiondiffusion (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(PCDit), 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_100Picardtrunc 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 .
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

1  4  7  57,120  63  36  30  30  \textcolorblack201  \textcolorblack29 
2  10  28  142,800  102  66  \colorblack59  54  \textcolorblack357  \textcolorblack48 
3  20  84  285,600  145  109  88  82  \textcolorblack553  \textcolorblack73 
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

1  4  10  57,120  26  22  12  11  \textcolorblack144  \textcolorblack13 
2  10  35  142,800  63  49  29  26  \textcolorblack300  \textcolorblack29 
3  20  84  285,600  145  109  88  82  \textcolorblack553  \textcolorblack73 
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

10  16  14  7  7  \textcolorblack186  \textcolorblack10 
20  36  27  17  16  \textcolorblack259  \textcolorblack17 
30  102  66  \textcolorblack59  54  \textcolorblack357  \textcolorblack48 
setup  0  1  2  3  6  

1  3  6  10  28  
0  12  21  43  63  
bGS  
CoV (%)  10  16  9  7  7  7 
20  36  19  15  17  17  
30  102  55  43  57  58  
ahGS  
CoV (%)  10  16  9  7  7  7 
20  36  19  14  16  16  
30  102  55  35  55  54 
Results for the first step of Newton iteration, which comes after six steps of Picard iteration, are summarized in Tables LABEL:tab:Re_100NewtonN–LABEL:tab:Re_100Newtontrunc. As above, the first three tables show the dependence on stochastic dimension (Table LABEL:tab:Re_100NewtonN), polynomial degree (Table LABEL:tab:Re_100NewtonP), and coefficient of variation (Table LABEL:tab:Re_100NewtonCoV). 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_100NewtonN and LABEL:tab:Re_100NewtonP). Nevertheless, checking results in Table LABEL:tab:Re_100Newtontrunc, 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.
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

1  4  7  57,120  73  42  32  32  \textcolorblack309  \textcolorblack40 
2  10  28  142,800  126  80  77  69  \textcolorblack546  \textcolorblack63 
3  20  84  285,600  235  170  128  151  \textcolorblack1011  \textcolorblack117 
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

1  4  10  57,120  27  23  12  11  \textcolorblack219  \textcolorblack32 
2  10  35  142,800  68  55  33  32  \textcolorblack450  \textcolorblack42 
3  20  84  285,600  235  170  128  151  \textcolorblack1011  \textcolorblack117 
MB  K  bGS  ahGS  \textcolorblackahGS(PCD)  \textcolorblackahGS(PCDit)  

10  17  15  8  8  \textcolorblack322  \textcolorblack28 
20  40  30  20  19  \textcolorblack379  \textcolorblack35 
30  126  80  77  69  \textcolorblack546  \textcolorblack63 
setup  0  1  2  3  6  

1  3  6  10  28  
0  12  21  43  63  
bGS  
CoV (%)  10  17  10  8  8  8 
20  40  22  18  20  20  
30  126  68  57  75  77  
ahGS  
CoV (%)  10  17  10  8  8  8 
20  40  22  17  19  19  
30  126  68  45  70  69 
Finally, we briefly discuss computational costs. For any preconditioner, each GMRES step entails a matrixvector product by the coefficient matrix. For viscosity given by a general probability distribution, this will typically involve a blockdense matrix, and, ignoring any overhead associated with increasing the number of GMRES steps, this will be the dominant cost per step. The meanbased preconditioner requires the action of the inverse of the blockdiagonal matrix . This has relatively small amount of overhead once the factors of are computed. The ahGS preconditioner without truncation effectively entails a matrixvector product by the block lowertriangular 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_100Picardtrunc and LABEL:tab:Re_100Newtontrunc, 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 meanbased 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 meanbased 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 blockdiagonal structure, and we do not discuss it here.
5 Conclusion
We studied the NavierStokes 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 GaussSeidel hierarchical preconditioning with approximation using the meanbased 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 meanvalue 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. SpringerVerlag, 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 lognormal 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. SpringerVerlag New York, Inc., New York, NY, USA, 1991. (Revised edition by Dover Publications, 2003).
 [11] L. Giraldi, A. Nouy, and G. Legrain. Lowrank approximate inverse for preconditioning tensorstructured linear systems. SIAM Journal on Scientific Computing, 36(4):A1850–A1870, 2014.
 [12] V. Girault and P.A. Raviart. Finite element methods for NavierStokes equations. SpringerVerlag, 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(89):607–616, 2000.
 [18] C. E. Powell and H C. Elman. Blockdiagonal preconditioning for spectral stochastic finiteelement systems. IMA Journal of Numerical Analysis, 29(2):350–375, 2009.
 [19] C. E. Powell and D. J. Silvester. Preconditioning steadystate 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 innerouter 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.