STOCHASTIC

. We study the time-dependent Navier–Stokes equations in the context of stochastic ﬁnite element discretizations. Speciﬁcally, we assume that the viscosity is a random ﬁeld given in the form of a generalized polynomial chaos expansion, and we use the stochastic Galerkin method to extend the methodology from [D. A. Kay et al., SIAM J. Sci. Comput. 32(1), pp. 111–128, 2010] into this framework. For the resulting stochastic problem, we explore the properties of the resulting stochastic solutions, and we also compare the results with that of Monte Carlo and stochastic collocation. Since the time-stepping scheme is fully implicit, we also propose strategies for eﬃcient solution of the stochastic Galerkin linear systems using a preconditioned Krylov subspace method. The eﬀectiveness of the stochastic Galerkin method is illustrated by numerical experiments.

1. Introduction.Models of mathematical physics are commonly based on partial differential equations (PDEs).In this study, we focus on the most popular PDE model in fluid mechanics, which is the Navier-Stokes equation [7,23].We consider a stochastic version of the model: we assume that the viscosity is given by a generalized polynomial chaos (gPC) expansion, we discretize the problem using spectral stochastic finite elements see, e.g., [13,23,25,39], and we wish to find the gPC expansion of the solution.The steady-state version of this problem was studied in [24,31,37], and our focus here is on the time-dependent counterpart.Our approach to time discretization is built on the fully implicit scheme with adaptive time-stepping strategy, which was developed for the deterministic Navier-Stokes equation by Kay et al. [21], see also [15].We extend their scheme in the stochastic Galerkin framework, and in particular we show that the physics inspired time-stepping strategy can be also adapted to this framework.The scheme is fully implicit, and so each time step entails a solve with the stochastic Galerkin matrix.This typically leads to very large systems of linear equations, for which use of direct solvers may be prohibitive, and therefore the method could potentially be quite computationally expensive.There are other approaches to time stepping see, e.g., [1,7,20] which may appear more appealing.Nevertheless, finally we also show that the iterative solvers by Sousedík and Elman [37], which are based on preconditioned Krylov subspace methods, are quite effective for the implicit time discretizations of the time-dependent Navier-Stokes problem as well.
Some aspects of the gPC methods for time-dependent problems were studied in literature see, e.g., [16,41,42].In particular, long-term integration was addressed by Gerritsma et al. [10], Heuveline, Schick and Song [19,34,36], Wilkins [38], Özen nad Bal [29,30], and most recently by Esquivel et al. [9], among others.Methods for flows exhibiting uncertain periodic dynamics were proposed, e.g., by Bonnaire et al. [2], Lacour et al. [22] and Schick et al. [33].These methods typically entail time-dependent or other variants of gPC expansions that are tailored to the changing character of the solution.Nevertheless, here we use a time-independent gPC basis, which turns out to be sufficient for the transient problems considered in our numerical experiments.Therefore, all these techniques can be viewed as complementary to the present study.We also note that Elman and Su [8] proposed a low-rank stochastic Galerkin solver based on monolithic (all-at-once) time discretization of the Navier-Stokes problem, however their scheme is based on a constant timestep.
Finally, we remark on possible interpretations of the Navier-Stokes problem with stochastic viscosity.In such case, the Reynolds number defined as where ν > 0 is the viscosity, U is the characteristic velocity and L is the characteristic length, is also stochastic.The possible interpretations of such setup are discussed by Powell and Silvester in [31]: for example, assuming fixed geometry, the stochastic viscosity is equivalent to Reynolds number being stochastic, which may correspond to a scenario when the volume of fluid moving into the channel is uncertain.The paper is organized as follows.In Section 2 we recall the algorithm for the deterministic problem, in Section 3 we formulate the algorithm for the stochastic problem using both the stochastic Galerkin and sampling methods, in Section 4 we report results of numerical experiments and provide details about the preconditioning of the Oseen problem, and finally in Section 5 we summarize and conclude our work.
2. Algorithm for the deterministic problem.We first recall the algorithm for the deterministic problem following Kay et al. [21].Let D ⊂ R 2 be a physical domain, and let T > 0 denote a stopping time.We wish to solve the time-dependent Navier-Stokes equation in D × [0, T ], where ( u, p) denote the fluid velocity and pressure, and ν ≡ ν(x) > 0 1 is the viscosity parameter, written as with boundary and initial conditions given on ) The initial velocity field is assumed to satisfy the incompressibility constraint, that is ∇ • u 0 = 0. We also assume that Γ N has nonzero measure so that the pressure is uniquely specified, and to this end we will use the outflow (do-nothing) boundary condition.We begin by recalling the implicit trapezoid rule (TR) as 1 The assumption ν(x) = const is the only difference from the setup in [21] in this section. where The nonlinear term is linearized as The linearization is based on extrapolation w n+1 − u n /k n+1 = u n − u n−1 /k n , from which we find Next, let (V D , Q D ) denote a pair of spaces satisfying the inf-sup condition and let V E be an extension of V D containing velocity vectors that satisfy the Dirichlet boundary conditions [4,7,14].The mixed variational formulation of (2.6)-(2.7)is: find for all ( v, q) ∈ V D × Q D .We note that p n+1 is not needed for subsequent time steps.Next, we recall the three ingredients of the algorithm as discussed in [21]: time integration, time-step selection and stabilization of the integrator.Time integration.Substituting u n+1 = u n + k n+1 d n into (2.9)-(2.10),rearranging and using D q (∇ • u n ) = 0, we get the so-called discrete Oseen problem: given u n , ∂ u n /∂t and the boundary update g := ( g n+1 − g n )/k n+1 , we first compute for all ( v, q) ∈ V D × Q D , and the TR velocity and acceleration are updated as

.13)
Time-step selection.The time step size is driven by the heuristic formula The local truncation error e n+1 is estimated by where the TR velocity u n+1 is compared with the AB2 velocity u n+1 * , which is computed using the explicit formula (2.16) There are three issues that need to be addressed: 1.The AB2 is not self-starting.To start the simulation we require a function u 0 with boundary data g 0 such that The initial acceleration (and pressure) is computed as follows: given the boundary update g := ( 12) is then constructed by setting n = 0 and defining ∂t , and its solution u 1 , p 1 is used to compute the acceleration at time t = k 1 as and allows to compute the AB2 velocity at the second time step.The startup is completed by switching on the time-step control at the third time step (k 1 = k 0 ). 2. Choice of initial time step.The strategy is to select a conservatively small value for k 0 , say 10 −8 .The time step then typically exhibits a rapid growth in the first few steps, roughly as k n+1 /k n = O (ε/eps) Stabilization of the integrator.The numerical stabilization is implemented using time-step averaging with the purpose to annihilate any contribution of the form (−1) n to the solution and its time derivative, which is invoked periodically every n * steps.For such a step the values of t * = t n and u * = u n are saved, we set t n = t n−1 + 1 2 k n , t n+1 = t * + 1 2 k n+1 and define the new "shifted" solution vectors as where d n is the TR update computed via (2.11)-(2.12).In our implementation, the parameter n * is fixed and the value is set to 10. Finite element formulation.We consider the discretization of the Oseen problem (2.11)-(2.12)by a div-stable mixed finite element method; in the numerical experiments we use Taylor-Hood elements see, e.g., [7].Let the bases for velocity and pressure spaces be denoted by {φ i } nu i=1 and {ϕ j } np j=1 , respectively.In matrix terminology, the Oseen problem at time step n entails solving a linear system where F n+1 is the velocity convection-diffusion matrix: a sum of the velocity mass matrix M, diffusion matrix A and convection matrix N n+1 , defined as where and w n+1 is computed from (2.8).The divergence matrix B is defined as (2.20) The right-hand side in (2.18) is constructed from the boundary data g n+1 , the computed velocity u n at the previous time level, and the acceleration ∂ u n ∂t .3. Algorithms for the stochastic problem.Let (Ω, F, P) represent a complete probability space, where Ω is the sample space, F is a σ-algebra on Ω and P is a probability measure.We assume that the randomness in the model is induced by a vector ξ : Ω → Γ ⊂ R m ξ of independent, identically distributed (i.i.d.) random variables ξ 1 (ω), . . ., ξ m ξ (ω), where ω ∈ Ω.Let B(Γ) denote the Borel σ-algebra on Γ induced by ξ, and µ denote the induced measure.The expected value of the product of measurable functions on Γ determines a Hilbert space where the symbol E denotes mathematical expectation.
In computations, we will use a finite-dimensional subspace T p ⊂ T Γ spanned by a set of multivariate polynomials {ψ (ξ)} that are orthonormal with respect to the density function µ, that is E [ψ k ψ ] = δ k , and ψ 1 = 1.This will be referred to as the gPC basis [40].The dimension of the space T p , depends on the polynomial degree.For polynomials of total degree p ξ , the dimension is n ξ = m ξ +p ξ p ξ .

Navier-Stokes equation with stochastic viscosity.
We use the same setup as in [37].Specifically, we consider that the expansion of viscosity is given as where ν (x) is a set of deterministic spatial functions, and index is related through a multi-index to the degrees of the random variables ξ 1 , . . ., ξ m ξ used in the construction of the gPC basis function ψ (ξ) see, e.g., [13,Section 2.4.3] or [39,Section 5.2].For simplicity, we will also assume that both the Dirichlet boundary conditions (2.3) and the initial condition (2.5) are deterministic.We seek a discrete approximation of the velocity in the form In literature it is sometime recommended to use a time-dependent gPC basis, that is ψ k (ξ, t), to keep the stochastic dimension low in long-time integration.However, this is a complementary strategy to the present study.Since it is not needed in our numerical experiments, we use only a time-independent gPC basis.

Stochastic Galerkin method.
The stochastic Galerkin formulation of problem (2.9)-(2.10)consists of using the expansion (3.2) and performing a Galerkin projection on the space T Γ using mathematical expectation in the sense of (3.1).That is, we seek velocity u n+1 ∈ T Γ ⊗ V E and pressure p n+1 ∈ T Γ ⊗ Q D for a given pair ( u n , p n ), such that and the stochastic counterpart of the discrete Oseen problem (2.11)-(2.13)is: given u n , ∂ u n /∂t and the boundary update and the TR velocity and the acceleration are updated as in (2.13).

Stochastic
Galerkin finite element formulation.The Galerkin projection leads to a large coupled system of equations with structure depending on the ordering of the unknown coefficientsts {u ik }, {p jk }.We will group velocity-pressure pairs for each k, the index of stochastic basis functions (and order equations in the same way), giving the ordered list of coefficients u 1:nu,1 , p 1:np,1 , u 1:nu,2 , p 1:np,2 , . . ., u 1:nu,n ξ , p 1:np,n ξ . (3.6) The discrete stochastic Oseen operator is built as follows.First, we set up the discrete components of the diffusion matrix using the expansion of viscosity (3.2) as Next, let w n+1 (x) denote the th term of the extrapolated velocity iterate (as in the expression on the right in (3.3) for k = ) at step n, and let Let n = max(n ν , n ξ ) and, if needed, define A = 0 for n ν < ≤ n and N n+1 = 0 for n ξ < ≤ n.Then in analogue to (2.19) define matrices which are incorporated into the block matrices These operators will be coupled with matrices arising from terms in T p , Combining the expressions from (3.10) and (3.11), using the ordering (3.6) yields the discrete stochastic Oseen system where ⊗ denotes the matrix Kronecker product.The entries of the vectors v and y are ordered as in (3.6).Note that H 1 is the identity matrix of order n ξ .Remark 3.2.With this ordering, which we used also in [37], the coefficient matrix contains a set of n ξ block 2 × 2 matrices of saddle-point structure along its block diagonal, given by This enables the use of existing deterministic solvers for the individual diagonal blocks.
We find it convenient to formulate the solvers in the so-called matricized format.To this end, we make use of isomorphism between R nxn ξ and R nx×n ξ determined by the operators vec and mat.Let n x = n u + n p and consider writing the solution of (3.12) using the ordering (3.6) . ., n ξ as in the expansions on the right in (3.3).Then we write v = vec(V), V = mat(v), where v ∈ R nxn ξ , V ∈ R nx×n ξ and the upper/lower case notation is assumed throughout the paper, so Y = mat(y), etc.Specifically, we define the matricized coefficients of the solution expansion where the column k contains the coefficients associated with the basis function ψ k .
In this setting, since (V ⊗ W) vec (X) = vec WXV T , the linear system (3.12) can be equivalently written as The time-step selection is driven by the formula (2.14), which we heuristically modify as follows.First, we run the deterministic solver with viscosity ν = ν 1 and record the set of time steps 0, t 1 , t 2 , , . . ., T .Then, we divide size of each interval [t n , t n+1 ] by n ξ , and we further round down the time-step size to the nearest power of 10.This procedure yields a sequence of time steps, which is then used for evolution of the stochastic Galerkin method.We note that an alternative strategy could be utilized by using the gPC coefficients corresponding to the mean velocity directly in formula (2.14), that is without an a priori run of the deterministic solver.

Sampling methods.
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 ξ (q) , which give realizations of the viscosity (3.2).That is, a realization of viscosity ν ξ (q) gives rise to deterministic functions u •, •, ξ (q) and p •, •, ξ (q) on D that satisfy the standard deterministic Navier-Stokes equations, and to corresponding finite-element approximations.
In the Monte Carlo method, the n M C sample points are generated randomly, following the distribution of the random variables ξ, and moments of the solution are obtained from ensemble averaging.In addition the coefficients in (3.3) could be determined at time t b using2 where for t b , b = 1, . . ., n b , we will consider an a priori set of time barriers, which is used in implementation to enforce all n M C instances of the deterministic solver to step through.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 [11,28].There are several ways to implement stochastic collocation to obtain the coefficients in (3.3).In the basic variant of the method, it is possible proceed either by constructing a Lagrange interpolating polynomial, or, in the so-called pseudospectral approach, by performing a discrete projection [39].We use the pseudospectral approach because it facilitates a direct comparison with the stochastic Galerkin method, and we refer, e.g., to [23] for an overview and discussion of integration rules.In particular, the coefficients in (3.3) are determined at time t b using a quadrature rule where ξ (q) and w (q) , q = 1, . . ., n q , are the collocation (quadrature) points and weights.Finally, we note that the other ways to perform stochastic collocation include the least-square approach and the compressed sensing approach see, e.g., [5,17,18,27].
4. Numerical experiments.We implemented the method in Matlab using the IFISS 3.5 package [6,35], and in this section we present results of numerical experiments for a model problem given by a flow around an obstacle.The geometry of the problem is shown in Figure 4.1.The discretization of the physical space consists of 12, 640 velocity and 1640 pressure degrees of freedom.The viscosity was taken to be a lognormal process, and its representation was computed from an underlying Gaussian random process using the transformation described in [12].That is, for = 1, . . ., n ν , ψ (ξ) is the product of m ξ univariate Hermite polynomials, and denoting the coefficients of the Karhunen-Loève expansion of the Gaussian process by g j (x) and η j = ξ j − g j , j = 1, . . ., m ξ , the coefficients in expansion (3.2) are computed as The covariance function of the Gaussian field, for points X i = (x i , y i ) ∈ D, i = 1, 2, was chosen to be where L x and L y are the correlation lengths of the random variables ξ i , i = 1, . . ., m ξ , in the x and y directions, respectively, and σ g is the standard deviation of the Gaussian random field.The correlation lengths were set to be equal to 25% of the width and height of the domain, i.e.L x = 3 and L y = 0.5.The coefficient of variation of the lognormal field, defined as CoV = σ ν /ν 1 where σ ν is the standard deviation, was set to either 1% or 10%.The stochastic dimension was m ξ = 2.The degree used for the polynomial expansion of the solution was p ξ = 3, and the degree used for the expansion of the lognormal process was 2p ξ , which ensures a complete representation of the process in the discrete problem [26].With these settings, n ξ = 10 and n ν = n = 28, and H is of order 10 in (3.12).For the mean value of viscosity we used ν 1 = 0.02, which corresponds to mean Reynolds number Re 1 = 100, and ν 1 = 6.67 × 10 −3 , which corresponds to mean Reynolds number Re 1 = 300.We note that the steadystate case was studied in [37], and the setup for the (deterministic) time-dependent problem is the same to [7, Chapter 10] except that the length of the channel was set to 12. Specifically, the initial condition for velocity was taken zero, and the Dirichlet boundary condition on the inflow ∂D Dir (the left side) was smoothly ramped up from zero to steady state as u (•, t) = 1 − e −5t w, where w is a Poiseuille (parabolic) flow profile, no-flow condition was prescribed on the top and bottom walls and natural 'donothing' condition was used on the outflow boundary (the right side).The initial time step was set to k 0 =10 −9 , and the problem was evolved from t 0 = 0 s to T = 10 s.The time-stepping method described in Section 2 was used for each sample of viscosity ν x, ξ (q) for both Monte Carlo and the stochastic collocation methods, and the method from Section 3.2 was used for the stochastic Galerkin method.In order to compare the gPC coefficients at the a priori chosen set of times, which we will refer to as time barriers, we prescribed the stochastic Galerkin solver to step through certain times, and we used the same set also for the Monte Carlo simulation in order to compare probability density function estimates of the velocity obtained by using all three methods.Specifically, we used time barriers t b = {0, 0.1, 0.2, 0.5, 1, 2, 5, 6, 8, 10}.The evolution of the time step for the two deterministic cases with mean Reynolds numbers Re 1 = 100 and Re 1 = 300 and for the stochastic Galerkin method is shown in Figure 4.2.We note that the heuristic used for the stochastic Galerkin methods yields the same time-step selection for both values of the Reynolds number.Specifically, only three steps with size 10 −9 are performed at the very beginning, then the step size increases to 10 −6 and eventually to 10 −5 for the most part of the first second.For the next two seconds it becomes 10 −4 and eventually 10 −3 for the rest of the time.
Now, let us consider first the case of Re 1 = 100 and CoV = 10%.Figure 4.3 shows the evolution of the gPC coefficients of the horizontal velocity, and the symbols and × represent the results of Monte Carlo and stochastic collocation at some of the time barriers.It can be seen that all methods are in agreement.Figure 4.4 shows the mean horizontal velocity, Figure 4.5 the variance of the horizontal velocity, and Figure 4.6 the variance of the vertical velocity at times 0.1s, 1s and 10s.From these figures it can be seen that the flow quickly evolves during the first second, and the later changes are relatively less dramatic.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 [7], and the variance of the horizontal velocity component is evolving to be concentrated in two "eddies" and it is larger than the variance of the vertical velocity component.In fact, it appears that all quantities are already at time 10s close to the steady state, see also Figures 4.13  The left panels show the pdf estimates of the velocity in x direction at points with coordinates (4.0100, −0.4339) (top), (4.0100, 0.4339) (bottom), where the variance of the velocity is relatively large cf. Figure 4.5.The right panels show the estimates at point (3.6436, 0) which is slightly downstream from the obstacle: the estimate in the x direction in the top panel and the estimate in the y direction in the bottom panel.The results were obtained using Matlab's ksdensity function.It can be seen that the changes of the mean values of the pdf estimates are relatively large during the first second, and then the uncertainty gradually increases and the supports of the pdf estimates grow as the solution evolves away from the deterministic initial condition and the effect of the stochastic viscosity becomes evident.Next, let us consider the case of Re 1 = 300 and CoV = 1%.Figure 4.8 shows the evolution of the gPC coefficients of the horizontal velocity.It can be seen that with increased Re 1 it takes more time for the flow to develop, including the stochastic components of the solution despite lower CoV than in the previous problem.Again, all methods are in agreement.vertical velocity, all at times 0.1s, 1s and 10s.The mean quantities are quite similar to what would be expected in the deterministic case, and the variances reflect on more complex behavior of the fluid at the higher value of Re 1 .Finally, Figure 4.12 displays evolution of the probability density function (pdf) estimates at the same set of points of the domain and times as Figure 4.7, and all three methods are again in agreement.Finally, we compare the results of the stochastic Galerkin method applied to the steady-state problem with mean Reynolds number Re 1 = 100 and CoV = 10%, which was studied by Sousedík and Elman in [37], and the results of the long-term integration at time 100s.Specifically, a comparison of the mean horizontal velocity is shown in Figure 4.13, and Figure 4.14 displays the variance of the horizontal velocity.By comparing the two figures, it can be seen that the results are virtually identical.(3.12) in each time step of the stochastic Galerkin method is a computationally expensive task.Therefore, use of a preconditioned Krylov subspace method may be preferred over a direct solver.To this end, we used the right-preconditioned flexible GMRES (fGMRES) method [32] with the so-called mean-based preconditioner M −1 1 : R −→ V, which entails solving a linear system where R and V are the matricized coefficients of the gPC expansions, cf.(3.13).Specifically, M −1 1 denotes an action of the pressure convection-diffusion (PCD) preconditioner, see [21,Section 3] and [7, Section 9.2.2], which is motivated by the block   where F n+1 1 is the matrix from (3.8), and X −1 is the pressure convection-diffusion term First, we used LU factorizations of the matrices from (4.2), which are updated in each time step.Since the solves with the (mean) matrix M 1 are thus exact, this illustrates the approximation properties of the mean-based preconditioner.Then, we also used the IFISS implementation of the PCD iterated preconditioner, in which the solves involving both F n+1 1 and A p = BT −1 B T , where T is the diagonal of the velocity mass matrix, are replaced by a single V-cycle of AMG using the IFISS default parameters, and the solve with the pressure matrix M p is effected by five Chebyshev iterations, see [7,Section 10.3].The construction of the matrix F n+1 p is described

Conclusion.
We studied the time-dependent Navier-Stokes equations with stochastic viscosity, which was given in terms of a polynomial chaos expansion.For this problem, we developed a stochastic Galerkin method with adaptive, mean-informed time stepping.We applied the method to a popular benchmark problem given by a flow around an obstacle, and we compared the solution of the time-dependent problem after the transient to that of the corresponding steady-state problem.Next, since the time-stepping scheme is fully implicit, a linear solve with the stochastic Galerkin matrix is required in each time step.Use of direct solvers may be prohibitive due to the large size of the systems, and in fact it is even not desirable to form the matrices explicitly.Therefore, we also formulated a preconditioner, which is used by the rightpreconditioned flexible GMRES method, and allows to solve the stochastic Galerkin systems efficiently.We studied two variants of the preconditioner.The first variant is based on exact factorization of the matrix corresponding to the underlying mean problem, and the second one was an iterated variant by means of an algebraic multigrid solver.In the numerical experiments we observed that the performance of the exact and iterated variants of the preconditioner was virtually identical, and only a couple of GMRES iterations were needed for convergence in all time steps.Therefore, the proposed stochastic Galerkin method is designed as a wrapper around an existing code for the corresponding deterministic problem, and in fact an efficient solver for the deterministic problem is the essential component also for the method presented in this study.Finally, we also compared the stochastic Galerkin solution with the stochastic collocation and Monte Carlo solutions, and we observed an excellent agreement for all problems studied in our numerical experiments.
and 4.14.A different perspective on the solution is given by Figure4.7, which displays evolution of the probability density function (pdf) estimates in several points of the domain at times 0.1s, 1s and 10s.

1 Fig. 4 . 2 .
Fig. 4.2.Evolution of the time-step size for the deterministic problems with Reynolds numbers Re 1 = 100 and Re 1 = 300 and for the stochastic Galerkin method (SG).

Fig. 4 . 3 .
Fig. 4.3.Evolution of the gPC coefficients corresponding to the horizontal velocity in terms of 2 -norm for mean Reynolds number Re 1 = 100 and CoV = 10%.The symbols and × represent the results of the Monte Carlo and stochastic collocation, respectively, at times 0.1s, 1s and 10s.

Figure 4 .
9 then shows the mean horizontal velocity, Figure 4.10 the variance of the horizontal velocity, and Figure 4.11 the variance of the

Fig. 4 . 8 .
Fig. 4.8.Evolution of the gPC coefficients corresponding to the horizontal velocity in terms of 2 -norm for mean Reynolds number Re 1 = 300 and CoV = 1%.The symbols and × represent the results of the Monte Carlo and stochastic collocation, respectively, at times 0.1s, 1s and 10s.

Fig. 4 . 14 .
Fig. 4.14.Variance of the horizontal velocity obtained using the stochastic Galerkin methods for the steady-state problem (top), and at time 100s (bottom) with mean Reynolds number Re 1 = 100 and CoV = 10%.