Elisa Schenone
email: elisa.schenone@inria.fr

Stéphane Veys
email: stephane.veys@imag.fr

Christophe Prud

HIGH PERFORMANCE COMPUTING FOR THE REDUCED BASIS METHOD. APPLICATION TO NATURAL CONVECTION

In this paper, we are interested in studying of the reduced basis methodology (RBM) applied to steady-state natural convection problems. The latter has applications in many engineering domains and being able to apply the RBM would allow to gain huge computation savings when querying the model the reduced model for many parameter evaluations. In this work, we focus on the order reduction of the model -in particular the handling of the non-linear terms, -as well as the design of the RBM computational framework and the requirements on high performance computing to treat 3D models using Feel++, a C++ open source library to solve partial differential equations. Numerical experiments are presented on a 2D and 3D models.

Introduction

Nowadays, in many application fields, engineering problems require accurate, reliable, and efficient evaluation of quantities of interest. Often, these quantities of interest depend on the solution of a parametrized partial differential equation where the -e.g. physical or geometrical -parameters are inputs of the model and the evaluation of quantities of interest are outputs -e.g. average values. -In a real-time or many-query context, the reduced basis method (RBM) offers a rapid and reliable evaluation of the input-output relationship (see [PRV + 02, VPP03, VPRP03, PP04, QRM11, RHP07] for the methodology) for a large class of problems.

In this paper, we are interested in studying the RBM applied to steady-state natural convection problems parametrized by the Grashof and Prandtl numbers. Natural convection has applications in many engineering domains and being able to apply the RBM would allow to gain huge computation savings when querying the model the reduced model for many parameter evaluations. In this work, we focus on the order reduction of the model -in particular the handling of the non-linear terms, -as well as the design of the RBM computational framework and the requirements on high performance computing (HPC) to treat 3D models. Even though the model considered remains simple with respect to industrial applications, we tackle some of the main difficulties namely order reduction and computational costs.

In order to solve Finite Elements (FE) or Reduced Basis (RB) problems, we use an open-source library called Feel++ for Finite Element Embedded Library and Language in C++ ( [PCD + 12, Pru06] ). Feel++ is a library to solve problems arising from partial differential equations (PDEs) with Galerkin methods, standard or generalized, continuous or discontinuous, from 1D to 3D, for low to high order approximations (including geometry). Among the many other Feel++ features, it provides a seamless programming environment with respect to parallel computing using MPI, see section 4.1. Feel++ enjoys an implementation of the RBM, see [START_REF] Daversin | A reduced basis framework : Application to nonlinear multiphysics problem[END_REF], which can deal with a wide range of problems: elliptic or parabolic models, coercive or noncoercive models, linear or non-linear models. It can handle coupled non-linear multiphysic problems such as the thermo-electric problems in [VDPT12, VCD + 12]. It is important that such an environment hides as many implementation details as possible and let the user worry only about his/her model and the high level aspects of the FEM and RBM.

The organisation of the paper is as follows: in section 1 we describe the 2D and 3D steady-state natural convection models; in section 2, we present the finite element discretization and the solution strategy; in section 3, we apply the RBM and focus in particular on the non-linear terms handling; in section 4, we present the computational framework for FEM and RBM as well as some implementation aspects; finally in section 5, we display some numerical experiments in 2D and 3D.

Problem setting

We start with the description of a standard natural convection model. We consider a heated fluid in a squared or cubical cavity, the fluid circulates towards the low temperature under the action of density and gravity differences. We introduce the adimensionalized steady-state incompressible Navier-Stokes equations coupled with the heat equation and we consider the problem in a two and a three dimensional tank, see e.g. Figure 1: find (u, p, T ) such that

                                       u • ∇u + ∇p - 1 √ Gr u = T e 2 , in Ω ∇ • u = 0 , in Ω u • ∇T - 1 √ GrP r T = 0 , in Ω u = 0 , on ∂Ω T = 0 , on Γ 1 ∂T ∂n = 0 , on ∂Ω \ (Γ 1 ∪ Γ 3 ) ∂T ∂n = 1 , on Γ 3 . (1)
where Ω ⊂ R d , d = 2, 3, u, p and T are respectively the adimensionalized velocity, pressure and temperature, Gr and P r are the Grashof and the Prandtl numbers, and e 2 is the inward-pointing normal vector of a border Γ 2 ⊂ ∂Ω. The 2D domain consists in a rectangular tank of height 1 and length W , in the 3D case we consider a rectangular cuboid of height 1, length W and depth 1. A heat flux is imposed on the "right" border Γ 3 while the temperature is fixed on the "left" border Γ 1 and the remaining walls are insulated. Similar boundary conditions apply in 3D. No-slip boundary conditions are set for the fluid velocity in the tank. The parameters are the Grashof and the Prandtl numbers and we consider the Nusselt number -the average temperature on Γ 3 -as the output. As the Grashof and/or Prandtl numbers increase the Nusselt is decreasing, see e.g. figure 5(a) or figure 6(a).

Finite element formulation

In this section, we write the weak Galerkin formulations associated to (1) and we propose an iterative method to solve this problem in case of FE space and a resolution with the RBM. The weak formulation associated to problem (1) writes

a(u, u, v) -b(p, v) + 1 √ Gr c(u, v) -d(T, v) = 0 , ∀v ∈ V b(q, u) = 0 , ∀q ∈ Q e(T, u, ξ) + 1 √ GrP r g(T, ξ) - 1 √ GrP r h(ξ) = 0 , ∀ξ ∈ Ξ (2) where V ≡ [H 1 0 (Ω)] d , Q ≡ L 2 (Ω), Ξ ≡ {ξ ∈ H 1 (Ω) s.t. ξ| Γ1 = 0}
, and we define the tri-linear forms a :

V × V × V → R and e : Ξ × V × Ξ → R as a(u, w, v) = Ω (w • ∇u) • v, ∀u, w, v ∈ V (3) e(T, v, ξ) = Ω (v • ∇T ) ξ, ∀v ∈ V, T, ξ ∈ Ξ (4) the bi-linear forms b : Q × V → R, c : V × V → R, d : Ξ × V → R and g : Ξ × Ξ → R as b(q, v) = Ω q∇ • v, ∀v ∈ V, q ∈ Q (5) c(w, v) = Ω ∇w : ∇v, ∀w, v ∈ V (6) d(ξ, v) = Ω ξ e 2 • v, ∀v ∈ V, ξ ∈ Ξ (7) g(T, ξ) = Ω ∇T • ∇ξ, ∀T, ξ ∈ Ξ (8)
and the linear operator h : Ξ → R as

h(ξ) = Γ3 ξ, ∀ξ ∈ Ξ. ( 9
)
We introduce now a FE discretization of (2). We define the discrete spaces V h ⊂ V , Q h ⊂ Q, Ξ h ⊂ Ξ and the associated Galerkin projection of the solution (u, p, T ) of (2) by

V h ≡ span{φ 1 , . . . , φ Nu } , u Nu j=1 u j φ j Q h ≡ span{ψ 1 , . . . , ψ Np } , p Np j=1 p j ψ j Ξ h ≡ span{ξ 1 , . . . , ξ N T } , T N T j=1 T j ξ j .
The discrete formulation of (2) now reads Nu i,j=1

u i u j a(φ i , φ j , φ k1 ) + Nu i=1 1 √ Gr u i c(φ i , φ k1 ) - Np i=1 p i b(ψ i , φ k1 ) - N T i=1 T i d(ξ i , φ k1 ) = 0 , k 1 = 1, . . . , N u Nu i=1 u i b(ψ k2 , φ i ) = 0 , k 2 = 1, . . . , N p (10) N T i=1 Nu j=1 T i u j e(ξ i , φ j , ξ k3 ) + 1 √ GrP r N T i=1 T i g(ξ i , ξ k3 ) - 1 √ GrP r h(ξ k3 ) = 0 , k 3 = 1, . . . , N T
Due to its strong non-linearities, when the Grashof or Prandtl numbers are high, a robust iterative method is required to solve this problem. We apply here a Newton Method. For a given parameter µ = (µ 1 , µ 2 ) = (Gr -1/2 , P r -1 ) and an initial guess (u 0 , p 0 , T 0 ), at each Newton sub-iteration n = 1, . . . , n max we look for

(u n+1 , p n+1 , T n+1 ) ∈ R Nu × R Np × R N T such that J(u n , p n , T n ; µ) (u n+1 , p n+1 , T n+1 ) -(u n , p n , T n ) = R(u n , p n , T n ; µ) (11)
where J = J(u, p, T ; µ) is the Jacobian matrix, R = R(u, p, T ; µ) is the residual vector.

In the case of problem (10) the terms of the Jacobian matrix can be easily calculated. For each row k 1 = 1, . . . , N u they write

J k1i (u, p, T ; µ) = Nu j=1 u j a(φ i , φ j , φ k1 ) + Nu j=1 u j a(φ j , φ i , φ k1 ) + µ 1 c(φ i , φ k1 ) , i = 1, . . . , N u J k1Nu+i (u, p, T ; µ) = -b(ψ i , φ k1 ) , i = 1, . . . , N p J k1Nu+Np+i (u, p, T ; µ) = -d(ξ i , φ k1 ) , i = 1, . . . , N T (12)
for all row

k 2 = N u + k, k = 1, . . . , N p J k2i (u, p, T ; µ) = b(ψ k , φ i ) , i = 1, . . . , N u J k2Nu+i (u, p, T ; µ) = 0 , i = 1, . . . , N p J k2Nu+Np+i (u, p, T ; µ) = 0 , i = 1, . . . , N T (13)
and for all line

k 3 = N u + N p + k, k = 1, . . . , N T J k3i (u, p, T ; µ) = N T j=1 T j e(ξ j , φ i , ξ k ) , i = 1, . . . , N u J k3Nu+i (u, p, T ; µ) = 0 , i = 1, . . . , N p J k3Nu+Np+i (u, p, T ; µ) = Nu j=1 u j e(ξ i , φ j , ξ k ) + µ 1 µ 2 g(ξ i , ξ k ) , i = 1, . . . , N T . (14)
In the same way, each term of the residual R(u, p, T ; µ) ∈ R Nu+Np+N T can be calculated by

R k (u, p, T ; µ) = -a(u, u, φ k ) + b(p, φ k ) -µ 1 c(u, φ k ) + d(T, φ k ) , k = 1, . . . , N u R Nu+k (u, p, T ; µ) = -b(ψ k , u) , k = 1, . . . , N p R Nu+Np+k (u, p, T ; µ) = µ 1 µ 2 h(ξ k ) -e(T, u, ξ k ) -µ 1 µ 2 g(T, ξ k ) , k = 1, . . . , N T Remark 2.1.
For high Gr and P r numbers the Newton method might be insufficient. We propose in that case to use the continuation algorithm 1. Note however that, for the range of Grashof and Prandtl, we didn't need to apply this continuation method.

Algorithm 1 Continuation strategy for high Gr and P r numbers.

Fix parameters Gr and P r Fix minimal values of parameters Gr min = 1 and P r min = 10 -2 Fix tolerance tol and max number of iteration n max for Newton Calculate number of intermediary parameters: N = max 1; max log(Gr/Gr min ) ), (log(P r/P r min ) for i=1:N do Logarithmic scale for intermediary parameters:

Gr(i) = exp{log(Gr min ) + i(log(Gr/Gr min ))/N ) P r(i) = exp{log(P r min ) + i(log(P r/P r min ))/N ) Fix Newton initial guess u 0 = u old while ||R|| ≥ tol or n ≤ n max do find u n s.t. J(u n-1 )(u n -u n-1 ) = R(u n-1 ) end while u old = u n end for return Solution u = u n

Reduced basis formulation

Let us now investigate a reduced basis formulation to solve the heat convection problem introduced above. We first propose a general approach that can be applied to any quadratic problem affine in parameters. The technique is based on the idea to store the more information as possible in the reduced space. So, tensors are introduced and projected as well as matrices and vectors. After an overview on this method we deal with application to the heat convection equations.

Reduced Basis for a general quadratic problem

We now turn to the generalization of the RBM to a quadratic problem, in abstract form it reads

a(u, u, v; µ) + b(u, v; µ) = f (v; µ), ∀v ∈ V (15)
where the solution is

u = u(µ) ∈ V , µ ∈ R Q indicates Q scalars parameters, V is an Hilbert space defined on a domain Ω ⊂ R d , with d = 1, 2, 3, a : V × V × V → R is a tri-linear form, b : V × V → R is a bi-linear form, and f : V → R is a linear operator.
We first introduce the FE approximation of (15) that leads to the RB approximation. The RB formulation of the general problem (15) applies then to (2).

Let us define a FE space V N ≡ span{v 1 , . . . , v N } ⊂ V , and the approximated solution of (15) as

u(µ) N i=1 u i (µ)v i
then the discrete formulation of (15) writes

N i,j=1 u i u j a(v i , v j , v k ; µ) + N i=1 u i b(v i , v k ; µ) = f (v k ; µ), ∀k = 1, . . . , N (16)
where the solution

u = u(µ) = [u 1 • • • u N ] T ∈ R N for each parameter µ ∈ R Q .
In order to describe the RB approximation applied to problems such as (15), we need to define some discrete objects. In particular, we introduce the tensor A = A(µ) ∈ R N ×N ×N , the matrix B = B(µ) ∈ R N ×N , and the

vector f = f (µ) ∈ R N defined by (A) ijk = a(v i , v j , v k ; µ), (B) ki = b(v i , v k ; µ), (f ) k = f k = f (v k ; µ)
for i, j, k = 1, . . . , N . Note that the tensor A can be considered as a vector of matrices, i.e. for each k = 1, . . . , N , we define the matrix

A k = A k (µ) ∈ R N ×N as (A k ) ij = (A) ijk , i, j = 1, . . .

.N . Using this notation (16) now reads

u T A k u + (Bu) k = f k , ∀k = 1, . . . , N (17)
where (Bu) k and f k are respectively the k-th term of vectors Bu and f . As described in Section 2 we can treat the non-linearity of the first term in (17) using a Newton Method. For a given parameter µ ∈ R Q , each iteration n = 1, 2, . . . of the Newton algorithm reads

J row(k) (u n ; µ)(u n+1 -u n ) = R k (u n ; µ), ∀k = 1, . . . , N (18)
where J row(k) (u n ; µ) ∈ R N and R k (u n ; µ) are respectively the k-th row of the Jacobian matrix

J = J(u; µ) ∈ R N ×N and the k-th term of the residual R = R(u n ; µ) ∈ R N defined by R k (u n ; µ) = f k -(u n ) T A k u n -(Bu n ) k , k = 1, . . . , N . (19)
For the sake of notation, the dependency of a, b, f , and all associated operators on the parameters µ is removed. Now, we just need to calculate the Jacobian matrix J = J(u; µ). We can write each of its elements as

J ki (u) = ∂a ∂u i (u, u, v k ) + ∂b ∂u i (u, v k ) - ∂f ∂u i (v k ) = = N l=1 N j=1 ∂ ∂u i u l u j a(v l , v j , v k ) + N j=1 ∂ ∂u i b(v j , v k )u j = = N l=1 N j=1 ∂ ∂u i u l A ljk u j + N j=1 ∂ ∂u i B kj u j = = N j=1 u j A jik + N j=1 A ijk u j + B ki . (20)
In general A is never assembled as it is of size N 3 , where N is the dimension of the underlying discretization space. If this is indeed the case for FE discretization, recall however that here we deal with reduced order approximation which enables the computation of A explicitly. Let us then introduce a reduced space V N = span{ϕ 1 , . . . , ϕ N }, with N N and define the projection

u = u(µ) ∈ R N of the solution u in V N as u = N i=1 u i ϕ i = Φ T u, (21) where
Φ = [ϕ 1 . . . ϕ N ] ∈ R N ×N , Φ ji = φi,j = (ϕ i , v j ) V N , the coefficients u i are u i = (u, ϕ i ) V N = (u, N j=1 φi,j v j ) V N = N j=1 u j φi,j and (•, •) V N is the scalar product associated to V N .
We observe that the Newton method defined in (18) is in fact generic with respect to the discrete spaces and that we can replace V N by V N which corresponds to projections of the terms in (18) onto V N . We now prove this statement.

Following the procedure of (21) we start by defining the projection of the source term

f as f = f (µ) ∈ R N f = Φ T f , (22)
then the matrix

B = B(µ) ∈ R N ×N as well is the projection of B into the reduced space V N B = Φ T BΦ, (23)
and finally we denote the reduced size tensor A = A(µ) ∈ R N ×N ×N . We now prove that it is in fact the projection of the tensor

A ∈ R N ×N ×N in V N A ijk = a(ϕ i , ϕ j , ϕ k ) = N l,m,h=1 φi,l φj,m φk,h a(v l , v m , v h ) = = N l,m,h=1 φi,l φj,m φk,h A lmh = N h=1 φk,h (Φ T A h Φ) ij = N h=1 φk,h ( A h ) ij (24)
i, j, k = 1, . . . , N , where

A k = Φ T A k Φ is the projection of A k ∈ R N ×N in V N for all k = 1, . . . , N . Let us define the tensorial product : R N × R N → R as a T (b T Ac) := N h=1 a h (b T A h c), ∀a, b, c ∈ R N , ∀A ∈ R N ×N ×N . (25)
Then

A ijk = N h=1 φk,h (ϕ T i A h ϕ j ) = ϕ T k (ϕ T i Aϕ j ), i, j, k = 1, . . . , N. (26)
Using equations ( 23) and ( 26) we obtain the projection of the Jacobian J ∈ R N ×N at each iteration n = 1, 2, . . ., and it reads

J ik ( u n ; µ) = ( A k | col(i) ) T u n + ( A k | row(i) ) u n + B ik ( 27
)
where

A k | col(i) = [ A jik ] N j=1 ∈ R N is the i-th column of the matrix A k for each k = 1, . . . , N and A k | row(i) = [ A ijk ] N
j=1 ∈ R N is the i-th row of the matrix A k for each k = 1, . . . , N . If we use the tensorial product defined by (25) we simply write

A k | col(i) = ϕ T k (ϕ T 1 Aϕ i ) • • • ϕ T k (ϕ T N Aϕ i ) A k | row(i) = ϕ T k (ϕ T i Aϕ 1 ) • • • ϕ T k (ϕ T i Aϕ N ) i, k = 1, . . . , N . Similarly, each term of the reduced residual R = R( u; µ) ∈ R N reads R k ( u n ; µ) = f k -( u n ) T A k u n -( B u n ) k , k = 1, . . . , N. (28)

Application to heat convection problem

Let us now apply this reduction technique to the natural convection problem introduced in Section 2. The FE formulation of each sub-iteration n = 1, . . . , n max of the Newton Method applied to equations (2) writes

J(u n , p n , T n ) (u n+1 , p n+1 , T n+1 ) -(u n , p n , T n ) = R(u n , p n , T n ) ( 29
)
with Jacobian and residual terms defined as in Section 2. In order to apply the technique described in the previous paragraph, we introduce matrices and tensors associated to (10). We define the tensors

A ∈ R Nu×Nu×Nu and E ∈ R N T ×Nu×N T as A = [A ijk , i, j, k = 1, . . . , N u ], (A k ) ij = A ijk = a(φ i , φ j , φ k ) E = [E ijk , i, k = 1, . . . , N T , j = 1, . . . , N u ], (E k ) ij = E ijk = e(ξ i , φ j , ξ k ) the matrices B ∈ R Np×Nu , C ∈ R Nu×Nu , D ∈ R Nu×N T , G ∈ R N T ×N T as B = [B ij , i = 1, . . . , N u , j = 1, . . . , N p ], B ij = b(ψ j , φ i ) C = [C ij , i, j = 1, . . . , N u ], C ij = c(φ j , φ i ) D = [D ij , i = 1, . . . , N u , j = 1, . . . , N T ], D ij = d(ξ j , φ i ) G = [G ij , i = 1, . . . , N T ], G ij = g(ξ j , ξ i )
and the vector

H ∈ R N T as H = [H i , i = 1, . . . , N T ], H i = h(ξ i ).
Then, using notation introduced in Section 3.1, the Jacobian matrix writes

J(u, p, T ; µ) =   J N l1 (u) + µ 1 C -B -D B T 0 0 J N l2 (T ) 0 J N l3 (u) + µ 1 µ 2 G   ( 30
)
where each element of the non-linear submatrices J N l1 (u), J N l2 (T ), J N l3 (u) are

J N l1 ki (u) = Nu j=1 ((A k ) ij + (A k ) ji )u j = ((A k | col(i) ) T + (A k | row(i) ))u , k, i = 1, . . . , N u J N l2 ki (T ) = N T j=1 (E k ) ij T j = (E k | row(i) )T , k = 1, . . . , N T , i = 1, . . . , N u (31)
J N l3 ki (u) = Nu j=1 (E k ) ji u j = (E k | col(i) ) T u , k, i = 1, . . . , N T
whereas each element of the residual vector of the Newton algorithm (29) writes

R k (u, p, T ; µ) = -u T A k u + B row(k) p -µ 1 C row(k) u + D row(k) T , k = 1, . . . , N u R k (u, p, T ; µ) = -(B col(k) ) T u , k = N u + 1, . . . , N u + N p (32) R k (u, p, T ; µ) =µ 1 µ 2 H k -u T E k T -µ 1 µ 2 G row(k) T
To apply efficiently the RB methodology, a key ingredient is the affine decomposition of the terms in the Newton Method which is readily available for our problem. The Jacobian matrix writes as

J(u, p, T ; µ) = Q J q=1 θ q J (µ)J q (u, p, T ) (33)
where for the considered example Q J = 4 and the coefficients θ J are

θ 1 J (µ) = µ 1 = 1 √ Gr , θ 2 J (µ) = µ 1 µ 2 = 1 √ GrP r , θ 3 J (µ) = θ 4 J (µ) = 1. ( 34
)
Each sub-matrix can be described using the notation introduced above

J 1 (u, p, T ) =   C 0 0 0 0 0 0 0 0   , J 2 (u, p, T ) =   0 0 0 0 0 0 0 0 G   , (35)
J 3 (u, p, T ) =   0 -B -D B T 0 0 0 0 0   , J 4 (u, p, T ) =   J N l1 (u) 0 0 0 0 0 J N l2 (T ) 0 J N l3 (u)   .
As to the residual, it is readily decomposed as

R(u, p, T ; µ) = Q R q=1 θ q R (µ)R q (u, p, T ) (36)
with Q R = 3 terms, we have

θ 1 R (µ) = µ 1 = 1 √ Gr , θ 2 R (µ) = µ 1 µ 2 = 1 √ GrP r , θ 3 R (µ) = 1 (37)
and

R 1 (u, p, T ) =   -Cu 0 0   , R 2 (u, p, T ) =   0 0 H -GT   , (38)
R 3 (u, p, T ) =   -(u T A k u) Nu k=1 + Bp + DT B T u -((u T E k T ) T ) N T k=1   .
The last step consists in projecting these equations in a reduced space. As shown for the general problem (15), the reduced Newton method at each iteration n = 1, . . . , n max writes

J( u n , p n , T n ) ( u n+1 , p n+1 , T n+1 ) -( u n , p n , T n ) = R( u n , p n , T n ) ( 39
)
where the ∼ indicates the projection onto the reduced space V N . The computation of vectors and matrices is done in a "classical" way, the terms are calculated in the FE space V N and projected into the reduced one by the matrix Φ defined by the basis function. Each basis function is the solution of problem (10) for a fixed parameter µ, possibly orthonormalized. For tensors A and E a different approach is needed because of their high dimensions.

Let us consider a general tensor A ∈ R N 3 defined by a trilinear form a : V 3 N → R, A ijk = a(v i , v j , v k ), i, j, k = 1, . . . , N as in problem (15). Its projection on a reduced space V N = span{ϕ 1 , . . . , ϕ N }, as shown in Section 3.1, is defined by A ∈ R N 3 , A ijk = ϕ T k (ϕ T i Aϕ j ), i, j, k = 1, . . . , N . We observe that A ijk = a(ϕ i , ϕ j , ϕ k ) where ϕ i , ϕ j , ϕ k ∈ V N .

Let us then define an hybrid tensor Λ ∈ R N ×N ×N whose elements are

Λ ijk = (Λ k ) ij = a(v i , v j , ϕ k ), (40)
i, j = 1, . . . , N , k = 1, . . . , N . We can then use this tensor to redefine the reduced tensor

A A ijk = a(ϕ i , ϕ j , ϕ k ) = N l,m=1 φi,l φj,m a(v l , v m , ϕ k ) = = N l,m=1 Φ li Φ mj (Λ k ) lm = (Φ T Λ k Φ) ij . (41)
So, we proved that

A k = Φ T Λ k Φ .
(42) That implies that we can calculate only N matrices Λ k , k = 1, . . . , N , and project them into the reduced space to obtain the reduced tensor A.

Computational framework

In this section we give an overview of the computational framework and its implementation to solve the heat convection problem introduced in section 2. All the routines described are part of the FE or RB frameworks of Feel++ library and available in Feel++ source code. We first introduce the main principles of Feel++, then we illustrate the FE implementation of problem (1) and finally we introduce the associated RB framework.

Feel++ principles and design

Feel++ is a C++ library that provides a clear and easy to use interface to solve complex PDE systems. It aims at bringing the scientific community a tool for the implementation of advanced numerical methods and high performance computing.

Feel++ relies on a so-called Domain Specific Embedded Language (DSEL) designed to closely match the Galerkin mathematical framework. In computer science DS(E)Ls are used to partition complexity. In Feel++ the DSEL splits low level mathematics and computer science on one side, and high level mathematics as well as physical applications to the other side. This difference between disciplines is reflected on users and developers tasks and allows easily improvements on both sides. Furthermore, it enables using Feel++ for teaching purposes, solving complex problems with multiple physics and scales or rapid prototyping of new methods, schemes or algorithms.

The DSEL on Feel++ provides access to powerful tools such as interpolation, with a simple and seamless interface, and allows clear translation of a wide range of variational formulations into the variational embedded language. Combined with this robust engine, it lies also arbitrary order finite elements, high order quadrature formulas and robust nodal configuration sets. The tools at the user's disposal grant the flexibility to implement numerical methods that cover a large combination of choices from meshes, function spaces or quadrature points using the same integrated language and control at each stage of the solution process of the numerical approximations.

In this paper, we use recent developments which allow to operate on large-scale parallel infrastructures. The general strategy used is parallel data framework using MPI and thanks to DSEL the MPI communications are seamless to the user: (i) we start with automatic mesh partitioning using Gmsh [GR09] (Chaco/Metis) -adding information about ghost cells with communication between neighbor partition;-(ii) the Feel++ parallel data structures such as meshes, (elements of) function spaces -create a parallel degrees of freedom table with local and global views; -(iii) and finally we use the library PETSc [BBB + 12b,BBB + 12a,BGMS97] which provides access to a Krylov subspace solvers(KSP) coupled with PETSc preconditioners such as Block-Jacobi, ASM, GASM. A complete description of the Feel++ high performance framework is available in the thesis [START_REF] Chabannes | Vers la simulation des écoulements sanguins[END_REF].

Remark 4.2. The last preconditioner is an additive variant of the Schwarz alternating method for the case of many subregions, see [START_REF] Smith | Domain decomposition: parallel multilevel methods for elliptic partial differential equations[END_REF]. For each sub-preconditioners (in the subdomains), PETSc allows to choose a wide range of sequential preconditioners such as LU, ILU, JACOBI, ML. Moreover, precondioner ASM or GASM can be used with or without an algebraic overlap. Other parallel preconditioners are available in PETSc but not used here. In particular we would like to mention the MUMPS direct parallel solver [START_REF] Amestoy | Multifrontal parallel distributed symmetric and unsymmetric solvers[END_REF]. We use it both as solver and preconditioner for iterative solves. FieldSplit preconditioners are also of notice for the applications we have: they allow to exploit the structure of block matrix.

Finite element model for the reduced basis framework

In this section we briefly introduce the organization of the RB framework of Feel++ (see figure 2), and then we see how the user deals with the FE model, needed to interface with the RB framework.

Reduced basis framework

The offline/online strategy developed in the RBM is implemented in the class CRBTrilinear. As the offline step of the method can be very expensive, scalar products resulting from the projection of matrices or vectors on the RB are saved in a database. To save objects in the database, we use the concept of "serialization" introduced by the set of libraries for the C++ programming language Boost. The class CRBApp is the driver for the RB framework. The class ParameterSpace manages the construction of the parameter space. First, if no existing database is found the offline step of RBM is run in order to build a basis, while if a database already exists the basis can be enriched or the algorithm passes to the online step.

CRBApp

To solve a linear or nonlinear problem with RB method in Feel++ the user has to define in his model the variational formulation of the problem. At the state of the art, the variation form must be defined in accord with the affine decomposition, future revisions of the code will introduce an automated affine decomposition. Note that using Feel++ the user model can be as well a 1D, 2D or 3D model. It is important to keep in mind that to interface with the RB framework the user has to provide only the FE model (i.e.parameter space, geometry, variational formulation), it corresponds to UserModel in figure 2 whose interface derives from CRBModelTrilinear.

The Feel++ RB framework support parallel architectures using the MPI technology. As in the FE framework an automatic mesh partitioning using Gmsh (Chaco/Metis) is computed, while every data associated to the reduced basis (scalars, vectors and dense matrices, parameter space samplings,...) are duplicated on each processor. However note that since the mesh is partitioned according the number of processors, finite element approximations and thus the reduced basis functions are in fact spread on all processors. Currently the basis functions are saved in the RB database with their associated partitioning. If they are required for visualization purposes or reduced basis space enrichment, the same data partition as in the initial computations must be used. Another particular attention must be paid to parameter space sampling generation: we must ensure that all processors hold the same samplings. To this end, there are generated in a sequential way by only one processor and then broadcasted to other processors.

Finite element model

Let us illustrate the implementation of the FE solution strategy for the problem (1) using Feel++, we give the main ideas of the solver function, with the continuation algorithm described above, and the computation of Newton method terms.

In listing 1, we display a snippet of code showing the code describing for a given µ, the solution process to retrieve (u(µ), p(µ), T (µ)): (i) compute the coefficients of the affine decomposition; (ii) assemble the linear terms; (iii) solve the non linear problem where updateJacobian and updateResidual are computing the jacobian and residual respectively during the Newton iterations. Note that updateJacobian need not be called at every iterations. Moreover the code is seamless with respect to geometrical dimension (2D or 3D) and parallel computing.

Listing 1. Implementation of algorithm 1 void solve ( parameter_type const & mu , element_ptrtype & T ) { this -> computeThetaQ ( mu ); // Update θ coefficients this -> update ( mu ); // Update affine decomposition of the linear terms

Flow profiles

First we look at the solutions obtained with RBM and compare them with the FEM ones. We investigate the relative error between the two resolutions for increased turbulent flow, i.e. high parameters.

Figure 3 shows 2D results for different Grashof values (1, 10 3 , 10 5 ). On the left there are RBM solutions, we show the velocity magnitude (top) and temperature (bottom) profiles. We observe faster flow for increased Grashof values as expected. On the right, relative errors for velocity (top) and temperature (bottom) are shown. As the flow is more turbulent for increasing parameters, the relative errors increase as well. Furthermore, errors is at least 2 × 10 -3 for velocity magnitude and 10 -2 for temperature profile, which are satisfactory values (we do not consider values for Gr = 1 because of it corresponds to a parameter used to build the basis, the error is as expected 10 -16 ).

We analyze the error in 2D and 3D simulations, respectively Figures 5(c) and 6(c), increasing Grashof values. We observe in both cases that the error increases for Grashof between 1 and 10 3 and it stabilizes for values greater than 10 3 . Also, we remark that the solution in a finite 3D tank for high Grashof value represents a flow with a complex pattern (Figure 4).

Performances

We now compare the RBM computational gain with respect to FEM. Figures 5(b) and 6(b) show the computational times for increased Grashof in both FEM and RBM cases, respectively in the 2D and 3D tanks. In both cases, we show the log-log plot of the computational times (in seconds) vs Gr.

As expected the FEM solution costs is several order more expensive than the RBM one. In particular, we have a factor 6 × 10 2 for small Grashof and we reach a factor 1.5 × 10 3 for high Grashof in 2D case, and a factor from 10 3 to 10 4 occurs in 3D case. Furthermore, we observe that the computational time increases in parameters in the case of FEM, while it is almost constant in RBM computations. In both 2D and 3D case we find few parameters that lead to a higher computational time in RB, this is due to increased number of online Newton iterations, however it still yield good results when inspecting the Nusselt number, see Figure 5(a) and Figure 6(a).

We remark that we do not apply the continuation algorithm neither for FEM nor for RBM. However it may become necessary for FEM resolution in case of non-convergence of the Newton solver, for higher parameter values or more complex geometry. As initial guess of the Newton Algorithm used to solve the online RB problems, the nearest known solution is used, whereas in FEM the initial guess is taken as the zero -although we could also use the nearest basis function as initial guess. -

Outputs

To conclude we look at the Nusselt number, i.e. the mean temperature on boundary Γ 3 = Ω ∩ {x = 1} (see Figure 1 for the 2D tank)

N u = Γ3 T dσ. (43)
This quantity decreases with increasing Grashof because of faster fluid flows that remove the heat for Γ 3 . In Figures 5(a) and 6(a), respectively 2D and 3D cases, the logarithmic output curves show that the RBM solutions follow the same behavior of the FEM ones with same N u values. These results confirm the good approximation obtained with RBM in both 2D and 3D simulations.

Conclusions and perspectives

We have presented a mathematical and computational treatment of the RBM applied to steady-state natural convection. This can be readily applied to other types of problems with quadratic non-linearities. The proposed method gives results that are accurate and efficient. In particular, the efficiency of this technique is remarkable in the case of flows with complex patterns both in 2D and 3D. From a framework point of view, the (offline) database handling raises interesting challenges when dealing with large parallel data. In terms of perspectives, a posteriori error estimation is a first step not only to assess the quality of the RB approximation but also to guide the RB space construction using greedy strategies [START_REF] Yano | A space-time petrov-galerkin certified reduced basis method : Application to the boussinesq equations[END_REF][START_REF] Veroy | A posteriori error bounds for reduced-basis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF]. More complex applications can be considered in particular including geometrical parameters. However it is foreseen that we may require hp-RBM approximations, see e.g. [START_REF] Eftang | A two-step certified reduced basis method[END_REF]. Finally, we dealt with the steady-state of the natural convection flow, the transient state is also of interest but requires much more involved mathematical and computational framework, see e.g. [START_REF] Yano | A space-time petrov-galerkin certified reduced basis method : Application to the boussinesq equations[END_REF].

Figure 1 .

 1 Figure 1. Geometry of the 2D model. Consider an extrusion of this geometry in 3D case is the extrusion in z axis of length 1.

Figure 2 .

 2 Figure 2. Class diagram for the Feel++ RB framework. Arrows represent instantiations of template classes.

  (a) Gr=1.e+00 (b) Gr=1.e+03 (c) Gr=1.e+05

Figure 3 .

 3 Figure 3. Comparisons of FEM and RBM solutions for Gr = (1, 10 3 , 10 5 ) and P r = 1. On the left velocity flow magnitude (top) and temperature profile (bottom) for RBM solutions. On the right velocity error (top) and temperature error (bottom).

  (a) mesh partitioning (24 processors) (b) temperature isosurfaces (c) Stream lines

Figure 4 .10Figure 5 .

 45 Figure 4. 3D computations for Gr = 1e7, P r = 0.1

Figure 6 .

 6 Figure 6. 3D computations for Gr ∈ [1; 1e6], P r = 1.

Acknowledgements

The authors wish to thank Vincent Chabannes, Ranine Tarabay and Céline Caldini-Queiros for fruitful discussions during the Cemracs 2012. Stéphane Veys and Christophe Prud'homme acknowledge the financial support of the project ANR HAMM ANR-2010-COSI-009. Finally the authors wish to thank the Cemracs 2012 and its organizers.

// nonlinear iterative solver , solve for U =(u, p, T ) backend () -> nlSolve ( _jacobian = updateJacobian , _residual = updateResidual , _solution = U ); } As mentioned earlier, at each Newton iteration, the Jacobian matrix and the residual vector are updated. We start with the linear terms which can be precomputed, see section 3.2. We show, as an example, the Jacobian sub-matrices assembly. The vector thetaAq[q] represents the coefficients θ J defined in (34), while the matrices Aq[0], Aq[1], Aq[2] are respectively the matrices J 1 , J 2 and J 3 defined in (35). The assembly of those matrices is displayed in listing 2.

Listing 2. Jacobian Linear terms assembly // Definition of mesh and FE space Xh = V h × Q h × Ξ h mesh_ptrtype mesh ; // N polynomial order Vh = Pch < N +1 , Vectorial >( mesh ); Qh = Pch <N >( mesh ); Xih = Pch < N +1 >( mesh ); Xh = Vh * Qh * Xix ; // Definition of functions element_type U ( Xh , " u " ); element_type V ( Xh , " v " ); // Velocity function and test function element_0_type u = U . element <0 >(); element_0_type v = V . element <0 >(); // Pression function and test function element_1_type p = U . element <1 >(); element_1_type q = V . element <1 >(); // Temperature function and test function element_2_type t = U . element <2 >(); element_2_type s = V . element <2 >(); // 1) Fluid equations // 1.1) Velocity diffusion :

integrate ( elements ( mesh ) , -idt ( t )*( trans ( vec ( cst (0.) , cst (1.0)))* id ( v ) ) ); // B . C . ... form2(Xh,Xh,M) builds a bilinear form X h × X h → R whose algebraic contribution will be stored in the matrix M. We remark that the finite element space Xh used to solve the problem is a composite space: Vh is a finite element space of degree 3 for vectorial functions (velocity), Qh is a scalar function space of degree 2 (pression) and Xih a scalar space of degree 3 (temperature). Again this is seamless for the user with respect to parallel computing and quite expressive with respect to the mathematical formulation.

Finally we add the contribution of the non-linear terms to the jacobian. The main operations required to this update are showed in listing 3. We remark that for a given µ the linear part of the jacobian stored in Jlin does not need to be updated, its implementation is computed once for all using the code displayed in listing 2.

Numerical Experiments

We present some numerical results: first we compare flow profiles obtained using FEM and RBM, and associated errors, in both 2D and 3D cases. Then, computational times and performances varying model parameters are show in both FEM and CRB cases are displayed. Finally, we compare the Nusselt number obtained in both cases.

Thanks to the Feel++ framework, the finite element and reduced basis models are available both in 2D and 3D.

Regarding the 2D case the FEM simulations refer to 6 × 10 4 degrees of freedom while in the 3D tank we have more than 2 × 10 5 dof. We consider polynomial basis functions of degree 3 for velocity and temperature variables, and degree 2 for pressure variable. Continuation algorithm (see Algorithm 1) is not used for FEM simulation and RBM as well, a simple Newton Method without parameters continuation is run. The RB used for all results presented here contains 28 elements non orthonormalized. The basis functions are solutions of the FE problem evaluated for Prandtl fixed to 1 while Grashof varies from 1 to 10 6 . Each 2D simulation is run in parallel on 10 processors, whereas each 3D simulation is run on 24 processors. In both cases we use the solver GMRES and the Additive Schwarz Method (GASM) as preconditioner.