Lorenzo Tamellini
email: lorenzo.tamellini@mail.polimi.it

Olivier Le Maitre

O Le Maître

♯ A Nouy

♭ Mox-Modellistica

Keywords: Uncertainty Quantification, Stochastic Navier-Stokes equations, Galerkin method, Model Reduction, Reduced Basis AMS Subject Classification: 60H35, 65C30, 76D05, 76M35, 60H30, 68U20

come

Model reduction based on Proper Generalized Decomposition for the Stochastic steady incompressible

Navier-Stokes equations

Model reduction based on Proper Generalized

Decomposition for the Stochastic steady incompressible Navier-Stokes equations *

Introduction

In recent years, functional approaches have been deeply investigated for the numerical solution of models driven by stochastic partial differential equations [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF]. These approaches consist in searching for a functional expansion of the random solution u on a basis of functions of a discrete set of random parameters modeling the input uncertainties. The solution u is thus a function defined on a parameter space Ξ equipped with a probability measure P , and with values in a certain function space V. Classical approximation methods consist in searching for an approximate M -terms expansion M k=1 u k Ψ k (ξ) of u(ξ), where the Ψ k are some suitable basis functions, typically polynomials or piecewise polynomials, and where the u k ∈ V are the coefficients that need to be computed. Approximate expansions can be computed using sampling-type approaches or Galerkin-type projection methods, these latter methods requiring the solution of a coupled system of M partial differential equations. For large-scale applications, the computation of these approximations becomes simply intractable.

In order to address this complexity, various model reduction methods have been proposed (see [START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF] for a short review). Model reduction methods based on nonlinear approximation aim at constructing an approximation of the parameterized solution u(ξ) under the form m i=1 u i λ i (ξ), where the u i and λ i constitute reduced bases of functions that are not fixed a priori but simultaneously determined using some suitable optimality criteria. These optimality criteria must be such that the m-term approximation is computable without any a priori information on the solution u.

A first class of model reduction methods, the so called "Reduced Basis" methods, define optimal approximations using a uniform norm on the parameter space [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF][START_REF] Maday | Global a priori convergence theory for reduced-basis approximation of single-parameter symmetric coercive elliptic partial differential equations[END_REF][START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]. For computational purposes, suboptimal approximations are introduced, using a greedy construction of deterministic approximation spaces. Reduced basis functions u i are progressively determined by the solution of successive deterministic problems associated with parameters values ξ i , i.e. u i = u(ξ i ), where a suitable error indicator is detected to be maximum. These methods have been applied to a large class of partial differential equations (see e.g. [START_REF] Veroy | Reduced-basis approximation of the viscous burgers equation: rigorous a posteriori error bounds[END_REF][START_REF] Quarteroni | Numerical solution of parametrized navierstokes equations by reduced basis methods[END_REF] for the application to Burgers and Navier-Stokes equations). Some convergence results have been recently obtained for a class of linear elliptic problems, under some regularity assumptions on the solution [START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF]. These approaches, initially introduced for parametric analyses, do not take into account the probability measure on the parameter space.

A second class of model reduction methods, known as Proper Generalized Decomposition methods (PGD), is based on the approximation of the weak solution of parametric/stochastic equations which is an element of a tensor product space V ⊗ S, where S is a space of functions defined on the weighted parameter space Ξ, typically S = L 2 (Ξ, P). It has been introduced in [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF] for the solution a class of linear stochastic PDEs, and then extended to other classes of stochastic PDEs (see e.g. [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | Generalized spectral decomposition method for stochastic non linear problems[END_REF][START_REF] Chevreuil | Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics[END_REF]). Different definitions of approximations have been proposed, which can be seen as generalized spectral decompositions (generalized Singular Value Decompositions). Based on the interpretation of the approximation problem as a nonlinear eigenproblem, several algorithms have been proposed, which are inspired from methods for the solution of eigenproblems, see [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF]. Greedy-type algorithms that construct the functions u i one after the other are of particular interest. Indeed, these algorithms only require the solution of successive deterministic problems. Note however that, unlike the aforementioned methods, these problems are not associated to a particular parameter value.

The PGD methods have also been successively applied to the solution of other high dimensional problems formulated in tensor spaces (see review [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF]). In particular, they have been used for the solution of high dimensional stochastic problems by further exploiting the tensor structure of stochastic function space [START_REF] Doostan | A least-squares approximation of partial differential equations with high-dimensional random inputs[END_REF][START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF]. General convergence results have been recently obtained for particular classes of elliptic problems [START_REF] Falco | A proper generalized decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF][START_REF] Cances | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF][START_REF] Falco | Proper Generalized Decomposition for nonlinear convex problems in tensor Banach spaces[END_REF]. Let us note that alternative solution strategies based on tensor approximation methods have also been proposed for the solution of high dimensional stochastic problems [START_REF] Ballani | A projection method to solve linear systems in tensor format[END_REF][START_REF] Khoromskij | Tensor-structured galerkin approximation of parametric and stochastic elliptic pdes[END_REF][START_REF] Matthies | Solving stochastic systems with low-rank tensor compression[END_REF]. These approaches are based on the use of classical tensor approximation methods within iterative solvers.

In this paper, we address the solution of the stochastic steady incompressible Navier-Stokes equations. Application of stochastic spectral methods to the Navier-Stokes equations, using Galerkin projection schemes, was first considered in [START_REF] Le Maître | A stochastic projection method for fluid flow. i. basic formulation[END_REF][START_REF] Le Maître | A stochastic projection method for fluid flow. ii. random process[END_REF][START_REF] Le Maître | Natural convection in a closed cavity under stochastic, non-Boussinesq conditions[END_REF][START_REF] Xiu | Modeling uncertainty in flow simulations via generalized Polynomial Chaos[END_REF], see also references in reviews [START_REF] Knio | Uncertainty propagation in cfd using polynomial chaos decomposition[END_REF][START_REF] Najm | Uncertainty quantification and polynomial chaos techniques in computational fluid dynamics[END_REF] and book [START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF]. Although successful, Galerkin methods for the stochastic Navier-Stokes equations are challenged by the dimension of the resulting non-linear problem and the need for adapted solvers [START_REF] Maître | A Newton method for the resolution of steady stochastic Navier-Stokes equations[END_REF]. Therefore, we propose in this work to apply the Proper Generalized Decomposition method to the stochastic Navier-Stokes equations. For this purpose, we extend to this nonlinear framework an algorithm that has been proposed in [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF] for the construction of the reduced basis of functions u i . This construction can be interpreted as an Arnoldi procedure for the solution of the associated nonlinear eigenproblem. Arnoldi iterations can be seen as a greedy procedure for the construction of a reduced approximation space. This algorithm has the remarkable property that for the construction of a m-dimensional reduced basis it only requires the solution of m deterministic PDEs that possess a classical structure, close to a deterministic incompress-ible Navier-Stokes problem. These deterministic problems can be handled by classical deterministic solvers, thus making the proposed algorithm a partially non intrusive method. The algorithm is applied to a divergence free formulation of the Navier-Stokes equations, yielding an approximation of the random velocity field on a reduced basis of divergence free deterministic velocity fields. A methodology is then proposed for the reconstruction of an approximation of the pressure field, the random velocity field being given. This approximation is defined through a minimal residual formulation of the Navier-Stokes equations. Two alternative methods are introduced for the construction of an approximation of the pressure. The first method is a direct application of a PGD algorithm to the minimal residual formulation of the Navier-Stokes equations, thus yielding to the construction of a convergent decomposition of the pressure. The second method, which is more computationally efficient, reuses as a reduced basis the deterministic pressure fields associated to the deterministic problems that were solved during the construction of the decomposition of the velocity field (i.e. the Lagrange multipliers associated with the divergence-free constraint).

The outline of the paper is as follows. In Section 2, the PGD method is presented in a general framework for the solution of parametric stochastic PDEs. In Section 3, we introduce the formulation of the steady incompressible Navier-Stokes equations and we detail the computational aspects of the application of the PGD. In Section 4, numerical examples illustrate the efficiency of the proposed method. Finally, the methodologies for pressure reconstruction are introduced in Section 5.

Proper Generalized Decomposition (PGD)

In this Section we introduce the weak-formulation of a generic problem with stochastic coefficients. We then shortly discussed the stochastic discretization using polynomial chaos expansion and the related Galerkin method. The Proper Generalized Decomposition method is then introduced and algorithms for its calculation are detailed.

Stochastic variational problem

Consider the following abstract deterministic variational problem:

Find u ∈ V such that a (u, v; φ) = b(v; φ), ∀v ∈ V, (1)
with V an appropriate vector space, φ the problem parameters, b(• ; φ) : V → R a linear form and a(•, • ; φ) : V × V → R a semi-linear form which is linear with respect to the second argument. The deterministic space V can be here either infinite or finite dimensional and is equipped with an inner product (•, •) V with associated norm • V . Note that if V has infinite dimension, it will have to be discretized at some point. However, to remain as general as possible, we delay the discussion on discretized spaces V to the next sections. In any case, we assume that problem (1) has a unique solution (depending on φ).

In this paper, we are interested in situations where the parameters φ of the problem are uncertain and therefore treated as random inputs. Let P := (Θ, Σ, µ) be an abstract probability space, where Θ is the set of random elementary events, Σ the σ-algebra of the events and µ a probability measure. For φ defined on P, we denote by φ(θ), θ ∈ Θ, a realization of the random parameters. The expectation of a generic random quantity h defined on P is denoted

E [h] := Θ h(θ) dµ(θ).
Let L 2 (Θ, µ) be the space of second-order real-valued random variables, equipped with the inner product (•, •) µ and associated norm

• L 2 (Θ,µ) , ∀(h, g) ∈ L 2 (Θ, µ), (h, g) µ := Θ h(θ)g(θ) dµ(θ), h L 2 (Θ,µ) = (h, h) 1/2 µ , so that h ∈ L 2 (Θ, µ) ⇔ h L 2 (Θ,µ) < +∞.
Since the parameters φ in equation ( 1) are random, the solution of (1), so denoted U , is also random and defined on P. It satisfies equation ( 1) almost surely, that is

Find U : Θ → V such that a.s. a (U (θ), v; φ(θ)) = b(v; φ(θ) ), ∀v ∈ V. ( 2
)
It will be further assumed that U ∈ V ⊗L 2 (Θ, µ), so that one can derive the fully weak variational form of the stochastic problem given by the following problem.

Stochastic problem.

Find

U ∈ V ⊗ L 2 (Θ, µ) such that A (U, V ; φ) = B(V ; φ), ∀ V ∈ V ⊗ L 2 (Θ, µ), (3)
with the forms A and B given by

A (U, V ; φ) := E [a (U, V ; φ)] = Θ a (U (θ), V (θ); φ(θ)) dµ(θ), B(V ; φ) := E [b(V ; φ)] = Θ b(V (θ); φ(θ)) dµ(θ).

Stochastic discretization

For computational purposes, numerical discretizations need to be introduced. These will concern both the deterministic space V, to be discussed in the following sections, and the stochastic space L 2 (Θ, µ), for which we rely on Polynomial Chaos (PC) expansions.

For the sake of simplicity, we restrict ourself to the case of PC approximations for a set of N independent identically distributed (i.i.d.) random variables, ξ = {ξ i , i = 1, . . . , N}, defined on P, with range Ξ and known probability density function dP(ξ). Any functional h : ξ ∈ Ξ → R is then a real-valued random variable and we have

E [h] = Θ h(ξ(θ)) dµ(θ) = Ξ h(y) dP(y).
In this context, we assume the knowledge of the random model parameters φ as a functional of ξ (see examples in the results sections), specifically

φ(θ) ≡ φ(ξ(θ)) a.s.
Since the model parameters are the only source of stochasticity in the problem, we have U (θ) ≡ U (ξ(θ)) for the solution of (2) with φ(ξ). In other words, the solution is computed in the probability space P(ξ) := (Ξ, Σ Ξ , dP), called the image space, instead of in the abstract space P. Further, we denote L 2 (Ξ, P) the space of second-order random variables, equipped with the inner product defined for (λ, β) ∈ L 2 (Ξ, P)

2 by

λ, β := Ξ λ(y)β(y) dP(y) = E [λβ] ,
and the associated norm

λ L 2 (Ξ,P) = λ, λ 1/2 = E λ 2 .
Next, we introduce an Hilbertian basis (complete orthonormal set) {Ψ 1 , Ψ 2 , . . . } of L 2 (Ξ, P), and denote by S M the subspace of L 2 (Ξ, P) spanned by the first M elements of the stochastic basis, that is

L 2 (Ξ, P) ⊃ S M := span {Ψ 1 , . . . , Ψ M }.
Any element λ ∈ L 2 (Ξ, P) can be approximated by λ M ∈ S M defined by the expansion

λ M (ξ) = M i=1 λ i Ψ i (ξ), lim M→∞ λ M -λ L 2 (Ξ,P) = 0.
Classically, the basis functions Ψ i are N-variate polynomials in ξ. Each standard measure P(ξ) over Ξ leads to a different classical polynomial family [START_REF] Xiu | The Wiener-Askey Polynomial Chaos for stochastic differential equations[END_REF], the case of ξ i standard Gaussian random variables corresponding to (normalized) Hermite polynomials [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF]. All developments below immediately extend to other types of stochastic basis, including piecewise polynomial approximations and hierarchical stochastic multi-wavelets. For spectral polynomial bases, a common truncature strategy is based on the maximal total degree of the basis functions retained in the construction of S M . Denoting No the maximal total degree, the dimension of

S M is dim(S M ) = M = (N + No)! N!No! ,
highlighting its combinatoric increase with both the number of random variables in ξ and the expansion degree No. Other possible construction strategies for S M have been investigated e.g. in [START_REF] Bäck | Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: a nu-merical comparison[END_REF].

Stochastic Galerkin formulation

The stochastic problem (3) can be recast in S M by means of Galerkin method, resulting in the following problem.

Discrete Stochastic Problem.

Find

U M ∈ V ⊗ S M such that A U M , V M ; φ = B(V M ; φ), ∀V M ∈ V ⊗ S M .
Inserting the PC expansion of the solution, U M = M i=1 u i Ψ i , in the previous equations results in a set of M coupled problems for the deterministic modes u i ∈ V of the solution [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF][START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF], namely

A M i=1 u i Ψ i , v l Ψ l ; φ = B(v l Ψ l ; φ), ∀v l ∈ V and l = 1, . . . , M. (4)
It is seen that dimension of the Galerkin problem is M times larger than the size of the original deterministic problem. Consequently, its resolution can be very costly, or even prohibitive, whenever N or No needs be large to obtain an accurate approximate U M of the exact stochastic solution. An additional difficulty appears when the form a is nonlinear in its first argument, making difficult the practical computation of the stochastic form A. These two difficulties call for improvement. First, regarding the dimensionality of the Galerkin problem, one can reduce complexity by relying on more appropriate expansion basis, e.g. by means of adaptive strategies and enrichment of polynomial basis (see e.g. [START_REF] Beck | On the optimal polynomial approximation of stochastic PDEs by Galerkin and collocation methods[END_REF][START_REF] Cohen | Analytic regularity and polynomial approximation of parametric and stochastic elliptic PDEs[END_REF][START_REF] Cohen | Convergence rates of best n-term Galerkin approximations for a class of elliptic sPDEs[END_REF]). However, adaptive approaches are complex to implement and often remains computationally intensive, while they do not address the difficulties related to nonlinearities. On the contrary, the PGD approaches discussed in the following aim at tackling the issues of dimensionality and, to some extent, are better suited to the reuse of deterministic code without special treatments of nonlinearities as a result. This latter point will be further discussed in the following.

PGD: principles

Let us go back to Formulation 3. The PGD method seeks for a separated representation of the solution U ∈ V ⊗ L 2 (Ξ, P) as

U (ξ) = i u i λ i (ξ),
where the u i ∈ V are the deterministic components and the λ i ∈ L 2 (Ξ, P) the stochastic components of the PGD. The m-terms PGD approximation of U , denoted U (m) , corresponds the truncated series

U (m) = m i=1 u i λ i ≈ U. ( 5
)
The objective is then to construct the expansion (5) to minimize the approximation error in some sense, without a priori selection of deterministic and stochastic components. PGD thus has to be contrasted with the classical Galerkin approach where the stochastic components, the Ψ i , are selected a priori, before the computation of the deterministic coefficients.

The simplest PGD algorithms determine the couples (u i , λ i ) ∈ V × L 2 (Ξ, P) one after the others. Specifically, assuming that U (m) has been already determined, let (u, λ) be the next couple of components. We here look for a correction uλ which lives in the manifold of rank-one elements in V ⊗ L 2 (Ξ, P). We here impose uλ to satisfy a Galerkin orthogonality with respect to the tangent manifold at uλ to the set of rank-one elements, which is defined by {uβ + vλ; β ∈ L 2 (Ξ, P), v ∈ V}. We therefore obtain the following necessary conditions for the definition of uλ:

Find (u, λ) ∈ V × L 2 (Ξ, P) such that A U (m) + uλ, uβ + vλ; φ = B(uβ + vλ; φ), ∀(v, β) ∈ V × L 2 (Ξ, P). (6)
For some classes of semilinear forms A, we can prove the existence of solutions uλ satisfying [START_REF] Cances | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF], see [START_REF] Falco | Proper Generalized Decomposition for nonlinear convex problems in tensor Banach spaces[END_REF]. Moreover, for some particular symmetric elliptic linear problems, the couples (u, λ) can be interpreted as left and right generalized singular vectors of U -U m , see [START_REF] Falco | A proper generalized decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF]. Among the solutions of ( 6), the best ones are selected by the algorithms described below that can be interpreted as algorithms for capturing approximations of the dominant singular vectors of U -U m . Note that for the present steady Navier-Stokes equations, the analysis of existence of solutions is still an open problem. Two coupled problems for u, λ can be derived from equation [START_REF] Cances | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF]:

Deterministic Problem. Find u ∈ V such that A U (m) + uλ, vλ; φ = B(vλ; φ), ∀v ∈ V. ( 7
)
Algorithm 1 Power method Solve deterministic problem: u ← D(λ; U )

1: U ← 0 [element 0 of V]
6: Normalize u: u ← u/ u V 7:
Solve stochastic problem: λ ← S(u; U )

8:
until (u, λ) converged 9:

U ← U + uλ 10: end for For λ given, we denote hereafter u = D(λ; U (m) ) the solution of deterministic problem [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF].

Stochastic Problem. Find λ ∈ L 2 (Ξ, P) such that A U (m) + uλ, uβ; φ = B(uβ; φ), ∀β ∈ L 2 (Ξ, P). (8)
Similarly, for u given, we denote λ = S(u; U (m) ) the solution of stochastic problem (8).

PGD: algorithms

The above interpretation of an optimal couple (u, λ) as a couple of dominant singular vectors of U -U (m) suggested to translate to the present situation techniques for the resolution of eigenvalues problems, like power-iteration or Arnoldi methods (see [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF]). Their application to scalar non linear problems has been thoroughly investigated in [START_REF] Nouy | Generalized spectral decomposition method for stochastic non linear problems[END_REF]. Note that these algorithms have also been investigated for other problems formulated in tensor product spaces, such as time-dependent partial differential equations [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF].

Power-Iterations

The power method for the computation of (u, λ) is stated in Algorithm 1. Note that the convergence criteria is not stated on the couple (u, λ) yielded by the power-type iterations is understood in a broad sense since u and λ may not converge individually (see [START_REF] Nouy | Generalized spectral decomposition method for solving stochastic finite element equations: invariant subspace problem and dedicated algorithms[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF] for discussion on the convergence of the iterations). In practice, only a limited number of iterations is performed. We also remark that λ and u have equivalent roles in the Algorithm, so that the normalization step at line 6 could be performed on λ rather then u.

The convergence of the resulting PGD obtained by the Power-Iteration algorithm can be improved by introducing an update of the stochastic components Solve deterministic problem: u ← D(λ; U )

8: Normalize u: u ← u/ u V 9:
Solve stochastic problem: λ ← S(u; U )

10: until (u, λ) converged 11:
Add u to its reduced basis:

W ← W ∪ {u} 12:
Add λ to its reduced basis: Λ ← Λ ∪ {λ} 13:

Solve update problem: Λ ← U(W )

14:

U ← l k=1 u k λ k 15: end for {λ 1 , . . . , λ m } after the determination of the m-th first couples. More specifically, given the deterministic components u 1 , u 2 , . . . , u m , the update problem consists in the solution of the following set of m coupled equations:

Update problem. Find λ 1 , . . . , λ m ∈ L 2 (Ξ, P) such that A m i=1 u i λ i , u l β; φ = B(u l β; φ), ∀β ∈ L 2 (Ξ, P), l = 1, . . . , m. ( 9
)
Denoting Λ (m) = {λ 1 . . . λ m }, the update problem is compactly written formally as

Λ (m) = U(W (m) ),
where

W (m) = {u 1 . . . u m } is called the reduced deterministic basis (of V).
The power-type algorithm with update is stated in Algorithm 2. Note that it is not necessary to solve the update problem (line 13 of Algorithm 2) at every step l. Moreover, it would be possible to update W instead of Λ. This would results in solving a Galerkin problem similar to the classical one, but with the stochastic basis {λ i } instead of the {Ψ i }

Arnoldi iterations

One disadvantage of Power-iterations-like methods is that they discard all the intermediate solutions within the repeat-until loops. The so-called Arnoldi algorithm is a possible solution to overcome such a "waste": the temporary solutions

l ← l + 1 8:
Solve deterministic problem u * ← D(λ; U )

9:
Orthogonalize

u * : u ← u * -l-1 k=1 (u k , u * ) V 10: if u V < ǫ then 11: l ← l -1 [stagnation of Arnoldi detected] 12:
Solve update problem: Λ ← U(W )

13:

U ← l k=1 u k λ k 14: else 15: Normalize u: u ← u/ u V 16:
Solve stochastic problem: λ ← S(u; U )

17:

Add u to its container:

W ← W ∪ {u} 18:
Add λ to its container: Λ ← Λ ∪ {λ}

19: if l = m then 20:
Solve update problem: Λ ← U(W )

21: U ← l k=1 u k λ k 22:
end if

23:

end if 24: end while are used to build a deterministic orthogonal basis W (m) , and then an update problem is solved to compute Λ (m) . The main advantage of this algorithm is therefore that it requires a lower number of resolutions for the determinstic and stochastic problems. The Arnoldi algorithm is stated in Algorithm 3.

Whenever the generation of deterministic modes stagnates into invariant subspaces (detected using the small positive parameter ǫ at line 10), an update step is performed. This update step can be interpreted as a deflation in the Arnoldi method. Note also that the update problems at lines 12 and 20 concern the whole stochastic components Λ generated so far, but one could as well perform a partial update considering only the Arnoldi subspace generated after the last detected stagnation.

Practical considerations

Obviously, also the algorithms above need a stochastic discretization. Again, we shall rely on PC expansions for the stochastic components and approximate the stochastic modes λ i in the finite dimensional S M by M k=1 λ k i Ψ k . Further, with this stochastic discretization, the stochastic problem [START_REF] Chevreuil | Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics[END_REF] and the update problem (9) translate into the Galerkin problems

A U (m) + u M k=1 λ k i Ψ k , uΨ l ; φ = B(uΨ l ; φ), l = 1, . . . , M, (10)
and

A m i=1 u i M k=1 λ k i Ψ k , u l Ψ j ; φ = B(u l Ψ j ; φ), l = 1, . . . , m and j = 1, . . . , M. (11)
For a given stochastic approximation space S M , one can expect the PGD solution U (m) to converge quickly to the Galerkin solution U M ∈ V ⊗ S M , with m ≪ M modes. This expectation comes from the fact that the PGD constructs the most relevant stochastic components λ i for the expansion, contrary to the Galerkin case where one chooses a priori the stochastic components (as the elements of the PC basis) and then seek for the solution in S M . Another point to be underlined in view of the above algorithms is that in each of them the computationally intensive steps are the resolution of the deterministic and, to a lower extent, the stochastic problems plus the update problems (optional in the Power-Iteration algorithm). As seen in ( 7) and ( 10) the size of the deterministic and stochastic problems are constant and equal to the dimension of the discretized spaces V and S M respectively; this is in general much lower than the size of the Galerkin problem which is the product of the two, with a significant complexity reduction as a result (provided that the number of systems to be solved is small enough). Concerning the update problem, we observe that its dimension is m×dim(S M ) so that if m is less than the dimension of the discretized space V the update problem is again much smaller in size than the Galerkin problem.

In addition, it will be shown in the following sections that for the Navier-Stokes equations the actual deterministic problems to be solved have structures very similar to the original Navier-Stokes equations, facilitating the re-use of existing deterministic codes, while implementing a Galerkin solver would require a greater implementation effort.

We also remark that instead of updating the stochastic components of the PGD solution, one could instead derive an update problem for the deterministic components {u i , i = 1, . . . , m}, which would in fact have the structure of the Galerkin problem in (4) but for the approximation in the stochastic space spanned by the {λ i } instead of the {Ψ i }. This alternative should be considered for problems where the dimension M of the stochastic space exceeds that of the discretized space V.

Navier-Stokes equations with uncertain parameters

We consider the bidimensional, steady, incompressible (constant density) Navier-Stokes equations on a bounded, simply connected domain Ω ⊂ R 2 with boundary ∂Ω. The dimensionless Navier-Stokes equations are

u • ∇u = -∇p + ∇ • σ(u) + f , ( 12a
)
∇ • u = 0, ( 12b
)
where u :

x ∈ Ω → R 2 is the velocity field, p : x ∈ Ω → R is the pressure field, f : x ∈ Ω → R 2
is the external force field and σ the viscous stress tensor. For a Newtonian fluid, σ in (12a) has for expression

σ(u) = ν 2 ∇u + ∇u T ,
where ν > 0 is the viscosity parameter (inverse of a Reynolds number), measuring relative influence of the inertial (nonlinear) and viscous (linear) contributions. Accounting for the mass conservation equation (12b), the Navier-Stokes equations reduce to

u • ∇u = -∇p + ν∇ 2 u + f , ( 13a
)
∇ • u = 0. ( 13b
)
These equations have to be complemented with boundary conditions; for simplicity, we shall restrict ourselves to the case of homogeneous Dirichlet velocity boundary conditions on ∂Ω,

u(x) = 0, x ∈ ∂Ω. ( 14
)
The case of non-homogeneous Dirichlet boundary conditions can be tackled by introducing a suitable affine space for the velocity, as shown in [START_REF] Nouy | Generalized spectral decomposition method for stochastic non linear problems[END_REF].

Functional framework

Next, we classically denote by L 2 (Ω) the space of functions that are square integrable over Ω. It is equipped with the following inner product and associated norm:

(p, q) := Ω pq dΩ, q L 2 (Ω) = (q, q) 1/2 .
We define the constrained space

L 2 0 (Ω) = q ∈ L 2 (Ω) : Ω q dΩ = 0 .
Then, let H 1 (Ω) be the Sobolev space of vector valued functions with all components and their first partial derivatives being square integrable over Ω, and H 1 0 (Ω) the constrained space of such vector functions vanishing on ∂Ω,

H 1 0 (Ω) = v ∈ H 1 (Ω), v = 0 on ∂Ω .
With the above notations, the Navier-Stokes system [START_REF] Falco | A proper generalized decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF] with boundary conditions [START_REF] Falco | Proper Generalized Decomposition for nonlinear convex problems in tensor Banach spaces[END_REF] then admits the following weak formulation.

Navier-Stokes equations.

Find (u, p) ∈ H 1 0 (Ω) × L 2 0 (Ω) such that c(u, u, v) + ν v(u, v) + d(p, v) = b(v), ∀v ∈ H 1 0 (Ω) ( 15
)
d(q, u) = 0, ∀q ∈ L 2 0 (Ω),
with the forms defined by

c(u, w, v) := Ω (u • ∇w) • v dΩ, v(u, v) := Ω ∇u : ∇v dΩ, d(p, v) := - Ω p∇ • v dΩ, b(v) := Ω f • v dΩ.
Pressure can also be formally suppressed in this weak formulation, by introducing the subspace of weakly divergence-free functions of H 1 0 (Ω), denoted hereafter H 1 0,div (Ω),

H 1 0,div (Ω) := v ∈ H 1 0 (Ω) : d(p, v) = 0, ∀p ∈ L 2 (Ω) .
Seeking u ∈ H 1 0,div (Ω), the weak form simplifies to the following problem. Divergence-free Navier-Stokes equations.

Find u ∈ H 1 0,div (Ω) such that c(u, u, v) + ν v(u, v) = b(v), ∀v ∈ H 1 0,div (Ω). (16)
Finally, we introduce the uncertain parameters. In this paper, we are concerned by situations where the external forcing f and viscous parameter ν are uncertain and, consistently with the previous sections, are seen as functions of a set of N random variables (e.g. normalized centered Gaussian random variables), ν = ν(ξ) and F = F (x, ξ). As a consequence, the divergence-free Navier-Stokes equation ( 16) has now a stochastic solution U (ξ). We can therefore state the following formulation:

Find U = U (ξ) : Ξ → H 1 0,div (Ω) such that c(U (ξ), U (ξ), V ) + ν(ξ) v(U (ξ), V ) = b(V ; F (ξ)), ∀ V ∈ H 1 0,div (
Find U ∈ H 1 0,div (Ω) ⊗ L 2 (Ξ, P) such that C(U , U , V ) + V ν (U , V ) = B(V ), ∀V ∈ H 1 0,div (Ω) ⊗ L 2 (Ξ, P). ( 17
)
The forms C, V ν and B are given by

C(U , W , V ) := E [c(U , W , V )] , V ν (U , V ) := E [ν v(U , V )] B(V ) := E [b(V ; F )] .
The previous formulation is ready to be discretized with the Stochastic Galerkin method, introducing the discretized stochastic space S M as in section 2.3. In practice, the divergence-free constraint is treated by adding a stochastic pressure field P (ξ), see e.g. [START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF]. Moreover, the size of the Galerkin problem is large, as all stochastic modes are coupled through the random viscosity and the non-linearity, so that efficient strategies for its resolution are needed, see for instance [START_REF] Maître | A Newton method for the resolution of steady stochastic Navier-Stokes equations[END_REF].

We will however base the following discussion on PGD on the formulation in H 1 0,div (Ω) ⊗ L 2 (Ξ, P) since we are looking for a PGD decomposition of U . We will return back to the issue of pressure later on.

PGD formulation

We now detail the deterministic, stochastic and update problems associated to the iterations of the PGD algorithms.

Deterministic problem

We here detail problem [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF]. We assume that a m-term reduced approximation U (m) = m i=1 u i λ i has been computed. For a given stochastic mode λ ∈ L 2 (Ξ, P), the associated deterministic mode u = D(λ; U (m) ) is defined by the following problem.

Find u ∈ H 1 0,div (Ω) such that C(λu, λu, λv)+C(λu, U (m) , λv) + C(U (m) , λu, λv) + V ν (λu, λv) = B(λv) -V ν (U (m) , λv) -C(U (m) , U (m) , λv), ∀v ∈ H 1 0,div (Ω).
For convenience and to stress the deterministic character of this problem we rewrite it as

Find u ∈ H 1 0,div (Ω) such that c(u, u, v) + c (u, v (m) c (λ), v) + c (v (m) c (λ), u, v) + ν v(u, v; λ) = b(v; U (m) , λ) ∀v ∈ H 1 0,div (Ω). ( 18
)
In the previous equation we have denoted

v (m) c (λ) = m i=1 E λ 2 λ i E [λ 3 ] u i , ν = E νλ 2 E [λ 3 ] b(v; U (m) , λ) = E [λ b(v; F )] E [λ 3 ] - m i=1 E [λνλ i ] E [λ 3 ] v(u i , v) - m i=1 m j=1 E [λλ i λ j ] E [λ 3 ] c(u i , u j , v) .
It is therefore seen that the structure of the deterministic PGD problem is essentially the same as the weak formulation of the deterministic incompressible Navier-Stokes equations, with a few remarkable differences. In particular: i) we have two new linear convective terms, associated with convective velocity v (m) c ; ii) the viscosity parameter is different, since its value is now ν = E νλ 2 / E λ 3 ; iii) the forcing term contains all the information about the previous modes which have been already computed. We further observe that we can always make ν > 0, by changing λ to -λ, owing to the homogeneity of the sought couple (λ, u).

As a result, the resolution of this problem can re-use existing deterministic flow solvers with minimal adaptations for the computation of the right-hand-side and the additional convection term. In addition, the enforcement of divergence free character of u can be achieved by introducing a deterministic Lagrange multiplier ∈ L 2 0 (Ω).

Stochastic problem

We now detail problem [START_REF] Chevreuil | Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics[END_REF]. Let us assume again that a m-term reduced approximation U (m) = m i=1 u i λ i has been computed. For a given deterministic mode u ∈ H 1 0,div (Ω), the associated stochastic mode λ = S(u; U (m) ) is solution of the following problem.

Find λ ∈ S M such that C(λu, λu, βu)+C(U (m) , λu, βu) + C(λu, U (m) , βu) + V ν (λu, βu) = B(βu) -C(U (m) , U (m) , βu) -V ν (U (m) , βu) ∀β ∈ S M
This is a quadratic equation for λ in weak form. We can highlight this by recasting the previous formulation as

Find λ ∈ S M such that E λ 2 β c(u, u, u) + m i=1 E [λλ i β] ( c(u i , u, u) + c(u i , u, u) ) + E [νλβ] v(u, u) = E [β b(u; F )] - m i,j=1 E [λ i λ j β] c(u i , u j , u) - m i=1 E [νλ i β] v(u i , u) ∀β ∈ S M . ( 19
)
To actually compute the PC expansion of λ in S M , λ = M k=0 λ k Ψ k , one has next to choose β = Ψ l in [START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF] and solve the following set of M quadratic equations in the coefficients

λ k : ∀ l = 1, . . . M , c(u, u, u) M k,k ′ =1 λ k λ k ′ E [Ψ k Ψ k ′ Ψ l ] + m i=1 ( c(u i , u, u) + c(u i , u, u) ) M k,k ′ =1 λ k λ i, k ′ E [Ψ k Ψ k ′ Ψ l ] + v(u, u) M k,k ′ =1 λ k ν k ′ E [Ψ k Ψ k ′ Ψ l ] = M k ′ =1 b ( f k ′ , u) E [Ψ k ′ Ψ l ] - m i,j=1 c(u i , u j , u) M k,k ′ =1 λ i, k λ j , k ′ E [Ψ k Ψ k ′ Ψ l ] - m i=1 v(u i , u) M k,k ′ =1 ν k λ i, k ′ E [Ψ k Ψ k ′ Ψ l ] ∀ l = 1, . . . M ,
where we have supposed that F admits a PC expansion,

F (x, ξ) = M k ′ =1 f k ′ (x) Ψ k ′ (ξ).

Update Problem

Finally, we detail the update problem [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF]. Given a m-term decomposition U (m) = m i=1 u i λ i , the update problem consists in recomputing all the m modes λ i by solving the following problem. [START_REF] Le Maître | A stochastic projection method for fluid flow. i. basic formulation[END_REF] In the present case, it consists in a system of m quadratic equations for λ i , all mutually coupled, but whose structure is close to the stochastic problem [START_REF] Le Maître | Spectral methods for uncertainty quantification[END_REF]. [START_REF] Le Maître | A stochastic projection method for fluid flow. i. basic formulation[END_REF], we end up with a system of quadratic equations for the coefficients λi,k , whose dimension is therefore m × M.

Find λ i ∈ S M , i = 1, . . . , m, such that C m i=1 u i λ i , m i=1 u i λ i , βu + V m i=1 u i λ i , βu = B(βu j ) ∀β ∈ S M , ∀j = 1, . . . , m
Denoting λ i = M k=1 λi,k Ψ k and taking β = ψ k , k = 1, . . . , M, in

Numerical results

In this Section we consider two test cases of increasing complexity and computational cost: in the first one the viscous parameter ν is the only uncertain parameter, while in the second one we consider both the viscous parameter and the forcing term as uncertainty sources. The aim of the tests is to compare the PGD approximation against the Galerkin solution, to assess the effectiveness of the method. All PGD approximations will be computed with the Arnoldi method described in Section 2.5.2.

As for the spatial discretization, we will consider a classical Spectral Element Method discretization, see e.g. [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF]. In particular, we will use a grid of Nu × Nu Gauss-Lobatto points for the approximation of the components of the velocity, while the pressure is approximated over a Nu -2×Nu -2 grid. The non linearity in the Navier-Stokes equation is solved with a preconditioned Quasi-Newton method, and at each step the linear system is solved with a GMRES solver. Once more we remark that the efficiency of the PGD method in determining the reduced approximation of U does not depend on the discretization method or Navier-Stokes solver considered, and any technique may be used.

Test 1: Random viscosity parameter

In the first test we consider a random viscosity ν given by

ν(θ) = ν ε + ν ′ (θ),
where ν ε > 0 and ν ′ (θ) has a Log-normal distribution with median value ν ′ > 0 and coefficient of variation C ν ′ ≥ 1. For these settings, the random viscosity can be expressed as

ν(θ) = ν ε + ν ′ exp (σξ(θ)) , σ := log C ν ′ 2.85 , ( 21
)
where ξ ∼ N (0, 1), ensuring that ν ′ ∈ [ν ′ /C ν ′ , ν ′ C ν ′ ] with a probability ≈ 0.995. Regarding the deterministic force field, it is well-known that force fields deriving from the gradient of a potential induce no flow for homogeneous boundary conditions. Therefore we consider the deterministic function ψ(x) and define f as

f = ∇ ∧ (0, 0, ψ) T , ( 22
)
so that ∇ ∧ f = (0, 0, -∇ 2 ψ) T . For simplicity, we restrict ourselves here to forcing terms having constant rotational Φ,

∇ ∧ f = (0, 0, Φ) T , ( 23
)
and a zero normal component on ∂Ω. This leads to the definition of ψ by

∇ 2 ψ = -Φ in Ω ψ = 0 on ∂Ω . ( 24
)
It is useful to further define the operator L : H -1 (Ω) → H 1 0 (Ω) that maps the forcing term Φ in [START_REF] Matthies | Solving stochastic systems with low-rank tensor compression[END_REF] to the corresponding solution, that is

L[Φ] = ψ . ( 25
)
The magnitude of the forcing term is fixed by Φ, which is hereafter set to Φ = 100 ν ′ to ensure that U Ω ≈ 1. The spatial structure of f is shown in Figure 1(a).

Galerkin solution

We start by setting ν ′ = 1/200, C ν ′ = 1.5, Nu = 51 and ν ε = 0.01ν ′ , and we consider the classical Galerkin Stochastic Projection method for the approximation of U . Guided by the expression of the viscosity in ( 21), we rely on a PC expansion of the solution using a single normalized Gaussian random variable ξ and corresponding Hermite PC basis. The Galerkin approximation is therefore sought as

U G (ξ) := No+1 k=1 u G k Ψ k (ξ), (26)
with No denoting the expansion order and Ψ k denoting the k-th degree Hermite polynomial in ξ. For this random viscosity distribution, a well converged solution is obtained for No = 10, as shown in the following discussion. The Galerkin solution for No = 10 is depicted in Figure 1(b)-1(d), showing the expected velocity field (that is the first mode of the Galerkin solution u G 0 , see Figure 1(b)), and the expectation and the standard deviation of the rotational of U G , see Figures 1(c) and 1(d). Plots in Figure 1 highlight the effect of nonlinearities. Indeed, since in the present situation the forcing term is deterministic and the viscosity parameter does not depend on x, if the nonlinear convective terms were neglected the solution of the resulting linear Stokes problem would be expressed as a product of a deterministic function times a stochastic factor, U (ξ) = α(ξ)u * . As a consequence mean and standard deviation of U would be equal to E(α)u * and to σ(α)u * respectively, and they would thus exhibit the same spatial structure. This is not the case here. Indeed, we observe in Figure 1(c)-1(d) that expectation and standard deviation field of the rotational of the velocity clearly exhibit different spatial patterns. In fact, the random viscosity has the strongest impact on the vorticity field along the boundary of the domain, where the shear stress is maximal and the uncertainty level reaches roughly 25%. Another stringent feature of the standard deviation of the vorticity field is the presence of detached structures along the boundary, that are created by the convective effects.

To better appreciate the complexity of the random flow field, as well as the converged character of the Galerkin solution for No = 10, the Karhunen-Loeve (SVD) decomposition of U G (ξ) is computed. Since the Galerkin solution is computed in a subspace S M , whose dimension is No + 1 = 11, its KL expansion is finite and writes as

U G (ξ) = No k=0 u G k Ψ k (ξ) = No+1 l=1 u G,KL l κ G l η l (ξ), κ G 1 ≥ κ G 2 ≥ • • • ≥ κ G No+1 ≥ 0, ( 27
)
where {u G,KL l } is an orthonormal set and E [η l η l ′ ] = δ ll ′ . Figure 2 shows the rotational of few KL modes u G,KL l : the plots show the increasing complexity with the mode index of the spatial structure of the rotational of the KL spatial modes. They also highlight the impact of the nonlinear convective term which induces a bending of these structures, due to the advection effects, which however possess the symmetries of the present problem.

Figure 2(b) shows the normalized spectrum, that is S l = κ G l / No+1 n=1 κ G n for l = 1, . . . , No + 1. It exhibits a fast decay, the 6-th normalized mode being 10 -5 times the first one, with essentially a uniform asymptotic decay rate except for the very last KL modes which are affected by the truncation of the stochastic basis.

PGD approximation

We next compute the PGD approximation of U , using the Arnoldi algorithm with ǫ = 0.01 and fixing the maximum rank of PGD to m = 15, and the KL decomposition of such PGD solution. We still use the same stochastic subspace S M as before. Figure 3 shows the expected velocity field (E[U (m) ]), and the expectation and standard deviation fields of the rotational of U (m) . The plots should be compared with those of the Galerkin solution shown in Figure 1, and the agreement is excellent.

The same conclusion arises when looking at the rotational of the KL spatial modes of the rank-15 PGD approximation, which are shown in Figure 4, and have to be compared with Figure 2. some of the spatial modes of the PGD approximation, and compares their rotationals with the rotationals of the corresponding KL spatial modes of the rank-15 PGD approximation. We observe that even if the PGD and KL modes exhibit similar structures, with finer details as the mode index increases, they are not the same. Indeed, the Arnoldi algorithm generates a sequence of orthogonal modes that are not the KL modes; however, the dominant modes of the KL expansion are well approximated in the successive spaces generated by the Arnoldi algorithm. Finally, we investigate the case where the viscosity parameter depends on more than one random variable. To do this, we modify the definition of ν from equation [START_REF] Le Maître | A stochastic projection method for fluid flow. ii. random process[END_REF] to

ν(θ) = ν ε + exp σ √ N ν Nν i=1 ξ i (θ) ,
with ξ i independent and normalized, centered, Gaussian random variables. This is clearly an over-parametrization of the problem, since indeed ξ

T (θ) = 1/ √ N ν Nν i=1 ξ i (θ)
is in turn a normalized, centered, Gaussian random variable; therefore ν truly has a unique stochastic dimension, such that the Navier-Stokes solution has the same intrinsic stochastic dimensionality ∀N ν ≥ 1.

It is found that the PGD approximation is quite insensitive to this overparametrization, thus proving to be able to capture the key features of the stochastic solution. This clearly appears in Figure 6, where we consider the rank-15 PGD approximations for problems with N ν = 1, 2, 3: here we compare the (normalized) norms of the PGD stochastic modes λ i (Figure 6(a)), and the (normalized) norms of the KL modes of the rank-15 PGD approximations (Figure 6(b)) for the three parameterizations tested. For the three values of N ν the decay of the λ i 's norm is essentially similar, although the deterministic ans stochastic problems are different and λ i are randomly initialized.s This is confirmed by the excellent agreement of the KL-spectra, which only differ for the last modes containing insignificant energy.

Test 2: Random forcing term

In the second test we consider also the forcing term as uncertain. To this end, we go back to equation [START_REF] Maître | A Newton method for the resolution of steady stochastic Navier-Stokes equations[END_REF] and take now Φ, the vertical component of the rotational of the force field, as a stationary Gaussian process with unit mean and standard deviation σ Φ > 0, characterized by the two point correlation function

C Φ (x, x ′ ) = E (Φ(x) -Φ 0 )(Φ(x ′ ) -Φ 0 ) = σ 2 Φ exp - x -x ′ L ,
where Φ 0 = 1 is the mean of Φ, L its correlation length, and xx ′ is the Euclidean norm in R 2 . The process admits the Karhunen-Loeve expansion

Φ(x, θ) = Φ 0 + ∞ i=1 Φ i (x)ξ i (θ),
where the ξ i are normalized uncorrelated Gaussian variables. Ordering the Karhunen-Loeve modes with decreasing norm Φ i L 2 (Ω) and truncating the expansion after the N f -th term results in the following approximation of the ex-

f 1 f 3 f 5 f 8
Figure 7: Spatial structure of some modes f i of the forcing term F for Test 2, as in equation [START_REF] Nouy | Recent developments in spectral stochastic methods for the numerical solution of stochastic partial differential equations[END_REF]. The length of the vectors have been scaled to improve for clarity.

ternal force field:

F (x, θ) ≈ F N f (x, θ) := f 0 + N f i=1 ξ i (θ) f i (x) , f i (x) = ∇ ∧    0 0 L[Φ i (x)]    . ( 28
)
We set L = 1, σ f / f 0 = 0.2, Nu = 35 and N f = 7. Figure 7 shows some of the modes f i of the forcing term. It is well known that as L decreases more and more KL modes are needed to represent accurately the forcing term. However, in this work we are not really concerned about the truncation error that stems from retaining only N f terms of the expansion, but only to show that the PGD method can handle such forcing terms in a natural way.

We consider again the viscous parameter ν as a lognormal random variable (N ν = 1), as in equation ( 21), and we set ν ′ = 1/100. This implies that the solution depends on N = N ν + N f = 8 random variables. The discrete probability space S M is selected setting No = 2, resulting in a set of M = 45 multivariate Hermite polynomials; within this setting, we compute the PGD solution up to rank-45, as well as the full Galerkin solution for validation purposes. In terms of computational cost, the rank-45 PGD solution requires the resolution of roughly 45 deterministic problems which amount to the core of the computational time.

Figure 8 shows the mean and standard deviation of the rotational of the PGD approximation for m = 45. We again observe the impact of the convective nonlinearities and the resulting high variability level, particularly pronounced along the domain boundary. The two first moments are in excellent agreement with the Galerkin solution (not shown).

Figure 9 shows some of the first PGD spatial modes. Contrary to the case of uncertain viscosity only, we know observe that the PGD modes u i have significant symmetry breaking, except for u 1 . This is again explained by the Arnoldi algorithm which aims at constructing orthogonal basis of dominant subspaces: the vectors spanning the subspaces don't necessarily reflect the symmetries of the solution in theses subspaces. However, the comparison of the KL spectra and modes of the PGD and Galerkin solutions (not shown) are in excellent agreement, proving that the successive Arnoldi subspaces effectively capture the dominant stochastic features of the solution.

Next, Figure 10 shows the decay of the total error norm with respect to the Galerkin solution. Specifically, we compute the normalized total error norm ǫ(m) through

ǫ 2 (m) := E U (m) -U G 2 Ω / E U G 2 Ω . ( 29
)
The error ǫ(m) is reported for ν ′ = 1/10, 1/50 and 1/100, the coefficient of variation being kept constant. As expected, when ν ′ decreases, the PGD rank increases to achieve a given error since the nonlinearity of the problem increases, with more and more complex stochastic features in the solution as a result. We also observe that the dimension of the successive Arnoldi subspaces tends to increase when ν ′ decreases, as shown by the separation between successive dots on the curves. However, for the lowest median viscosity value the rank-45 PGD solution has a relative error less than 10 -5 , and is obtained at a fraction of the Galerkin computational cost. In addition, we also monitor the convergence of the relative norm of the last λ i added during the Arnoldi procedure, which can be considered as a very naive error estimator. Indeed, the deterministic modes u i are normalized, and therefore whenever λ i is small the correction λ i u i becomes negligible. Such an error estimator, although very rough, turns out to be quite effective; we will reconsider it in the next Section. Finally, we repeat the convergence analysis for the PGD-Galerkin error on a second test-case, in which we consider N f = 14 terms in the KL expansion of f , thus obtaining a problem with N = 15 random variables. We consider again a polynomial stochastic space of order No = 2, whose dimension is M = 861. For such a problem a full Galerkin approach may be very computationally demanding both in terms of CPU-time (a few hours) and memory occupation, therefore we only monitor the convergence of the error indicator proposed above. The results are very encouraging, since we obtain errors of at least 10 -6 using only 60 modes. We remark that this roughly corresponds to the computational cost of the resolution of about 60 Navier-Stokes problems and a few sets of coupled quadratic equations for the gPCE coefficients of λ i .

Residual computation and pressure reconstruction

At this point, it is crucial to devise an error estimator to stop the PGD procedure as soon as the reduced solution is close enough to the exact solution in H 1 0,div (Ω) ⊗ S M . The most natural approach would be a stopping criterion involving the evaluation of the norm of the residual of the Stochastic Navier-Stokes equation ( 17) associated to the m-terms reduced solution U (m) in the discretized space H 1 0,div (Ω) ⊗ S M . The Arnoldi algorithm would then be stopped as soon as such residual becomes lower than a given tolerance in a suitable norm.

In practice, computing the residual of the Navier-Stokes equations in their divergence-free formulation [START_REF] Knio | Uncertainty propagation in cfd using polynomial chaos decomposition[END_REF] is not a convenient operation. Therefore, we go back to the weak deterministic Navier-Stokes equations [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF] and introduce the Stochastic Velocity-Pressure Navier-Stokes equations.

Find U ∈ H 1 0 (Ω) ⊗ S M , P ∈ L 2 0 (Ω) ⊗ S M such that C(U , U , V ) + V ν (U , V ) + D(P , V ) = B(V ) ∀V ∈ H 1 0 (Ω) ⊗ S M , ( 30
) D(Q, U ) = 0 ∀Q ∈ L 2 0 (Ω) ⊗ S M ,
where D(Q, V ) is defined as the expected value of the bilinear form d(•, •) appearing in [START_REF] Ghanem | Stochastic Finite Elements: A Spectral Approach[END_REF],

D(Q, V ) = E [d(Q, V )] .
Computing the residual for the velocity-pressure formulation is an affordable task, but at this point the PGD algorithm has not provided us with an approximation of the stochastic pressure yet. Hence, we now introduce a procedure to recover the pressure P (m) associated to the m-terms PGD solution U (m) .

Computing such approximation will introduce some computational overhead, but one could be interested in an approximation of the pressure anyway. We stress that the notation P (m) does not refer to an m-terms approximation of P , but to a generic approximation of P given the m-terms reduced approximation of U .

Pressure computation

For easiness of presentation, let us define

N (W , V ) := C(W , W , V )+V ν (W , V )-B(V ), ∀V , W ∈ H 1 0 (Ω)⊗S M , ( 31
)
and let V , W denote the scalar product in H 1 0 (Ω)⊗S M . Inserting the m-terms PGD velocity U (m) and the corresponding pressure P (m) into the Stochastic Velocity-Pressure Navier-Stokes equations [START_REF] Nouy | Proper Generalized Decompositions and separated representations for the numerical solution of high dimensional stochastic problems[END_REF] we have

N (U (m) , V ) + D(P (m) , V ) = R (m) , V ∀V ∈ H 1 0 (Ω) ⊗ S M , (32a) D(Q, U (m) ) = 0 ∀Q ∈ L 2 0 (Ω) ⊗ S M , (32b)
where R (m) denotes the residual of the momentum equation (32a), R (m) ∈ H 1 0 (Ω) ⊗ S M . Note that the continuity equation (32b) has no residual; indeed, all the deterministic modes in U (m) are divergence-free, being solutions of the deterministic problem [START_REF] Le Maître | Natural convection in a closed cavity under stochastic, non-Boussinesq conditions[END_REF]. Equation [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF] states that the residual R (m) is a function of the pressure P (m) . Hence, we propose here to define P (m) as the minimizer of R (m) in some prescribed norm. To be more computationally oriented, we next derive the problem for P (m) in the discrete case.

Let us denote by V h ⊂ H 1 0 (Ω) the finite dimensional velocity space, and with Π h ⊂ L 2 0 (Ω) the finite dimensional pressure space. Upon the introduction of the bases for V h and Π h defined in [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF] and that will be used in the results sections, we can identify any element W h ∈ V h ⊗ S M with the coordinates in the respective basis W h ∈ R dim(V h ) ⊗ S M , and similarly any element

Q h ∈ Π h ⊗ S M with Q h ∈ R dim(Π h ) ⊗ S M ;
in other words, W h (ξ) and Q h (ξ) are vectors whose components are functions of ξ, belonging to the subspace S M ⊂ L 2 (Ξ, P). Equation (32a) can therefore be recast as a semidiscrete equation in R

dim(V h ) ⊗ S M , N (m) h (ξ) + D T P (m) h (ξ) = R (m) h (ξ), (33)
with

N (m) h (ξ), R h (ξ) ∈ R dim(V h ) ⊗ S M , P (m) h (ξ) ∈ R dim(Π h ) ⊗ S M
, and D ∈ R dim(Π h )×dim(V h ) the deterministic discrete divergence operator. Next we define the residual norm as

R (m) h (ξ) 2 = R (m) h (ξ) 2 R dim(V h ) ⊗S M = E R (m) h (ξ) 2 R dim(V h ) ,
Thus, using [START_REF] Quarteroni | Numerical solution of parametrized navierstokes equations by reduced basis methods[END_REF], we obtain that the pressure minimizing

R (m) h (ξ) is the so- lution of D D T P (m) h (ξ) = -D N (m) h (ξ) , (34)
Note that D D T is a deterministic operator, and equation ( 34) is well-posed if V h and Π h verify the inf-sup condition. Moreover, computing the PC expansion of P (m) h (ξ), that is

P (m) h (ξ) = M k=1 P (m) h,k Ψ k (ξ), with P (m) h,k ∈ R dim(Π h ) and Ψ k (ξ) ∈ S M Hermite polynomials, results in a set of M uncoupled problems D D T P (m) h,k = -D N (m) h,k .

Note that N (m)

h,k has to be computed, using the projection N

(m) h,k = E [ N (m) h (ξ)Ψ k (ξ)],
since the stochastic vector N (m) h (ξ) derives from a non-linear combination of the PGD solution, hence its PC expansion is not immediately available.

Even if we can take advantage of this by factorizing the operator DD T only once to improve the computational efficiency (e.g. with a LU, ILU or Cholesky factorization), the overall cost may be demanding if the discrete stochastic space S M is large: indeed, it would require the resolution of M independent systems. One could then apply a PGD procedure to obtain an approximation of the stochastic pressure

P (m) h (ξ) = m ′ k=1 P (m) h,k γ k (ξ), (35)
with

P (m) h,k ∈ R dim(Π h
) and γ k (ξ) ∈ S M generic functions, using any of the Algorithms illustrated in Section 2.5 to solve [START_REF] Rozza | Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations[END_REF]. Note that the PGD approximation of P may in general use m ′ = m modes. Further savings can be achieved by using as deterministic modes for the pressure the Lagrange multipliers obtained during the resolution of the deterministic steps of the PGD decomposition of U (m) ; note that in this case m = m ′ . h (ξ) , the residual being computed using the pressure reconstructed with a PGD approach; iv) "λ norm" denotes the normalized norm of λ i , that is λ i / i λ i 2 .

Numerical results

In the previous section we have proposed two ways of computing an approximation of the pressure field: a "fully reduced approach" in which we use a PGD method to compute both the deterministic and the stochastic modes of the decomposition [START_REF] Veroy | Reduced-basis approximation of the viscous burgers equation: rigorous a posteriori error bounds[END_REF], and a "partly reduced approach" in which the deterministic modes P (m) h,k of ( 35) are taken to be the Lagrange multipliers resulting from the solution of the deterministic problems during the Arnoldi iterations. In both cases, the obtained pressure approximation will be then used to compute the residual R (m) h (ξ) through equation [START_REF] Quarteroni | Numerical solution of parametrized navierstokes equations by reduced basis methods[END_REF], and the norm R (m) h (ξ) will be used as a stopping criterion for the Arnoldi method.

We now aim at assessing the performances of these two stopping criteria, along with the one proposed in Section 4, i.e. the monitoring of λ i . Such criterion may be reasonable whenever one is not at all interested in pressure reconstruction, or willing to reconstruct P (m) h,k only once, after a satisfying approximation W (m) h,k has been computed. The convergence of the proposed quantities for Test 2 is shown in Figure 11. The residual computed by recycling the Lagrange Multipliers is slightly worse than the one computed after having reconstructed the pressure with a PGD approach. We observe that residual norms clearly overestimate the error in solution by 1-2 orders of magnitude, hence representing a quite restrictive criterion for the convergence of the method. On the other hand, the norms of the λ i appear closer to the true error, but slightly underestimating it, hence representing an "optimistic" criterion.

Conclusions

In this work we have investigated the resolution with a PGD method of the steady-state Navier-Stokes equations with uncertain forcing term and Reynolds number. Based on the interpretation of PGD as a model reduction method associated with a generalized spectral decomposition, iterative algorithms have been proposed for the progressive construction of reduced bases of approximation. Different algorithms have been presented, which are inspired from solution methods for solving eigenproblems. In particular, we have employed an Arnoldi method in our numerical simulations. The proposed strategy relies on the fact that the Navier-Stokes problem can indeed be recast as a problem for the velocity field only (in the space of divergence-free functions), and the pressure can be reconstructed in a second step. Future works should therefore investigate the application of PGD methods to situations in which such "problem reduction" is not possible.

A key feature of such PGD methods is that the computation of the deterministic and stochastic modes of the solution is decoupled: this allows remarkable savings with respect to the standard Galerkin technique, both in terms of computational complexity and coding effort required. We have indeed shown that because of such decoupling it is possible to reuse any existing solver with minimal adaptations for the computation of the deterministic modes. In the case of the Navier-Stokes equations, these adaptations simply entail the modification of the convective velocity, the viscosity parameter and the forcing term. The stochastic and update problems can also be solved with available software, since they amount to systems of quadratic equations.

The convergence of the PGD approximation of the velocity to the full Galerkin solution has been investigated in different numerical settings. In all the considered cases, the PGD is able to provide reasonable approximations of the full Galerkin solution with a limited number of modes, (10 -4 at least with approximately 20 modes in the test we have considered, see figure 10), thus with a smaller computational cost compared to the solution of the full Galerkin problem.

However, care has to be taken in the reconstruction of a reduced pressure: the mathematical formulation of this problem is non-trivial, and we have addressed this topic only for the discrete problem, proposing different approaches with different computational costs and achievable accuracy. This is certainly an aspect worth a deeper investigation in future works, as well as the extension of the proposed technique to non-steady problems.

2: for l in 1 , 2 ,

 12

Algorithm 2 4 :

 24 Power method with update 1: U ← 0 2: W ← {} [initialization of the reduced basis in V] 3: Λ ← {} [initialization of the reduced basis in L 2 (Ξ, P)] for l in 1, 2, . . . , m do

  Ω), for a.e. ξ ∈ Ξ, whose fully weak counterpart can be written immediately as Stochastic Navier-Stokes problem.

Figure 1 :

 1 Figure 1: Reference Galerkin approximation for Test 1.

Figure 2 :

 2 Figure 2: KL expansion of the reference Galerkin approximation of Test 1.

Figure 4 (

 4 b) shows the matching between the spectra of the two KL decompositions, again showing good agreement between the solutions. Figure5shows

Figure 3 :

 3 Figure 3: Rank-15 PGD approximation of Test 1.

Figure 4 :

 4 Figure 4: KL decomposition of the rank-15 PGD approximation for Test 1.

Figure 5 :

 5 Figure 5: Test 1. PGD modes number 1, 6, 11 (left), their rotational (center) and the rotational of the corresponding KL modes of the rank-15 PGD approximation (right).

  spectrum, N=1 PGD-KL spectrum, N=2 PGD-KL spectrum, N=3 (b) normalized norms of KL modes

Figure 6 :

 6 Figure 6: Comparison of the rank-15 PGD approximations of Test 1 with N ν = 1, 2, 3.

Figure 8 :Figure 9 :

 89 Figure 8: Rank-45 PGD approximation for Test 2. mode 1 mode 5 mode 8 mode 10

  , err. est. ν=1/100, ||u G -u PGD || ν=1/100, est. , err. est. ν=1/50, err. est. ν=1/100, err. est.

Figure 10 :

 10 Figure 10: Test 2: convergence of error ǫ(m) and of the relative norm of λ m with respect to the rank (m) of the PGD approximation for ν ′ = 1/10, 1/50, 1/100, considering N = 8 (left) and N = 15 (right) random variables.

Figure 11 :

 11 Figure11: convergence of the quantities proposed as stopping criterion for the PGD method with respect to the number of modes m. i) "error" denotes the normalized PGD-Galerkin errorU (m) -U G / U G ; ii) "LM-residual" denotes the normalized norm of residual R (m) h (ξ) / R (0)h (ξ) , the residual being computed using the Lagrange Multipliers as deterministic modes for the pressure; iii) "PGD-residual" denotes the normalized norm of residual R (m)h (ξ) / R (0)

Le Maître and A. Nouy are partially supported by GNR MoMaS (ANDRA, BRGM, CEA, EdF, IRSN, PACEN-CNRS) and by the French National Research Agency (Grants ANR-08-JCJC-0022 and ANR-2010-BLAN-0904). L. Tamellini is partially supported by the Italian grant FIRB-IDEAS (Project n. RBID08223Z) "Advanced numerical techniques for uncertainty quantification in engineering and life science problems".