Cécile Daversin
email: cecile.daversin@lncmi.cnrs.fr

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

Christophe Trophime
email: christophe.trophime@lncmi.cnrs.fr

Christophe Prud

A REDUCED BASIS FRAMEWORK: APPLICATION TO LARGE SCALE NON-LINEAR MULTI-PHYSICS PROBLEMS

In this paper we present applications of the reduced basis method (RBM) to large-scale non-linear multi-physics problems. We first describe the mathematical framework in place and in particular the Empirical Interpolation Method (EIM) to recover an affine decomposition and then we propose an implementation using the open-source library Feel++ which provides both the reduced basis and finite element layers. Large scale numerical examples are shown and are connected to real industrial applications arising from the High Field Resistive Magnets development at the

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 -e.g. average values -are outputs. 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 [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF],Veroy et al., 2003a[START_REF] Veroy | A posteriori error bounds for reducedbasis approximation of parametrized noncoercive and nonlinear elliptic partial differential equations[END_REF][START_REF] Prud | Reduced-basis output bounds for approximately parameterized elliptic coercive partial differential equations[END_REF],Quarteroni et al., 2011, Rozza et al., 2007] for the methodology) for a large class of problems.

In this paper, we are interested in studying the RBM applied to large scale non-linear multi-physics parametrized partial differential equations requiring not only a robust mathematical framework but also a HPC-enabled computational framework. We propose an implementation of the reduced basis method and the extensions to non-linear and non-affinely parametrized problems. Other implementations are available, such as [A.T. Patera andG.Rozza, 2007, D.J.Knezevic and[START_REF] Knezevic | A high-performance parallel implementation of the certified reduced basis method[END_REF]. The reduced basis methodology is suited to develop efficient strategies to tackle design and optimization in industrial context however, to our knowledge, it has not (yet) been used effectively in this context. Typical activities of a design department require the ability to efficiently perform parametric studies and sensitivity analysis to improve and guide engineers in their daily work. In particular, safety is an essential ingredient in industrial investigations which involves to take into account all eventual uncertainties on input parameters. The reduced basis framework provides a valuable tool for design only if it may be applied seamlessly from small to large scale applications and may take care of non-linear quantities. We will show how Feel++ framework -which includes these features -may be used to help engineers in their design process. Numerical examples are given in the context of the High Fields Resistives magnets development at the french Laboratoire National des Champs Magnétiques Intenses. To our knowledge these examples are among the first to show applications of the reduced basis methodology on industrial problems that lead to actual realizations.

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++ ( [Prud'Homme et al., 2012a, Prud'homme, 2006] ). 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 2.1. Feel++ enjoys an implementation of the RBM which can deal with a wide range of problems: elliptic or parabolic models, coercive or non-coercive models, linear or non-linear models. As mentioned earlier, 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 RB mathematical framework on which we build the RB computational framework described in section 2. In section 3, we present numerical experiments on some of the non-linear multi-physics applications of LNCMI.

Mathematical framework for reduced basis

This section describes the mathematical framework of the reduced basis method (RBM) and ingredients we need. We first introduce the outline of RBM on elliptic linear problems with affine dependence in parameters. The non-linearity that comes accross the problems we focus on can lead to a non-affine dependence on input parameters. The Empirical Interpolation Method (EIM) is a good way to manage this by recovering affine dependence in parameters. Furthermore the rapid evaluation of the output offered by RBM, the framework gives a reliability guarantee for the results. This uses Successive Constraints Method (SCM) as an ingredient to compute efficiently the lower bound of the coercivity constant of a bilinear form, in order to have an a posteriori error estimation. Finally, the application of RBM on elliptic non-linear problems by means of the tools previously introduced is more precisely described.

1.1. Elliptic linear problems with affine dependence in parameters 1.1.1. Preliminaries

Let Ω be a suitably regular bounded spatial domain in R d (for d=1,..,3). Denoting L 2 (Ω) the space of square integrable functions over Ω, we have H 1 (Ω) = {u|u ∈ L 2 (Ω), ∇u ∈ L 2 (Ω) d } and H 1 0 (Ω) = {u ∈ H 1 (Ω)|u ∂Ω = 0}. From an Hilbert space X ≡ H 1 0 (Ω) ν -or more generally H 1 0 (Ω) ν ⊂ X ⊂ H 1 (Ω) ν , where ν = 1 (respectively d) for a scalar (respectively vector) field, we define X N a finite element approximation space of (typically very large) dimension N .

General problem settings

Let u(µ) be the solution of a parametrized Partial Differential Equation (PDE) with respect to the input parameter p-vector µ ∈ D, where D ⊂ R p is the parameter space. We are interested in the evaluation of an output of interest s(µ) ∈ R which can be expressed as a functional of a field variable u(µ) :

s(µ) = (u(µ)) , (1)
for a suitable linear operator (.). The variational formulation of the PDE consists in finding u

(µ) ∈ X(Ω) such that a(u(µ), v; µ) = f (v; µ) ∀v ∈ X(Ω) , (2)
where a(., .; µ) and f (., µ) are respectively bilinear and linear forms associated with the PDE. An important ingredient of the RBM is the development of an efficient offline/online strategy. To this end, a and f and must depend affinely in µ that is to say that there exists positive integers Q a , Q f and Q such that a(., .; µ), f (.; µ) and (.; µ) can be expressed as

                       a(u, v; µ) = Qa q=1 θ q a (µ)a q (u, v) ∀u, v ∈ X(Ω), ∀µ ∈ D , f (v; µ) = Q f q=1 θ q f (µ)f q (v) ∀v ∈ X(Ω), ∀µ ∈ D , (v; µ) = Q q=1 θ q (µ) q (v) ∀v ∈ X(Ω), ∀µ ∈ D , (3)
where 3 and 1.4 propose solutions when the affine decomposition are not readily available.

θ q a : D → R , 1 ≤ q ≤ Q a , θ q f : D → R , 1 ≤ q ≤ Q f and θ q : D → R , 1 ≤ q ≤ Q are µ-dependent functions. Sections 1.

Reduced basis method

We now turn to the construction of the reduced basis approximation. For given a µ ∈ D, we start with the finite element (FE) discretization of problem (1)-( 2) which consists in evaluating

s N (u N (µ)) = (u N (µ); µ) , (4)
where

u N (µ) ∈ X N satisfies a(u N (µ), v; µ) = f (v; µ) ∀v ∈ X N .
(5) For a given positive integer N max we introduce a nested sequence of reduced basis approximation spaces W Npr , 1 ≤ N ≤ N max . Note that W Npr is a N-dimensional subspace of X N , and that N max is very small in comparison to N . Let S N = {µ 1 , ..., µ N }, for each µ i ∈ S N we define the set of solutions of (5)

S u N = {u N (µ i ) , ∀µ i ∈ S N } .
(6) The application of the Gram-Schmidt process with respect to the (., .) X inner product to elements of the set S u N gives mutually (., .) X -orthonormal basis functions ξ pr n , 1 ≤ n ≤ N . The reduced basis space W Npr is then defined as

W Npr = span{ξ pr n , 1 ≤ n ≤ N } . (7)
The reduced basis solution u N (µ) can be expressed as

u N (µ) = N j=1 u N j (µ) ξ pr j . (8)
Now we use the affine parameter dependence to construct an efficient offline/online strategy. Choosing test functions as

v = ξ pr i , i = 1, . . . N , u N (µ) then satisfies N j=1 Qa q=1 θ q a (µ)a q ξ pr j , ξ pr i u Nj (µ) =   Q f q=1 θ q f (µ)f q (ξ pr i )   , (9)
which can also be written in a matrix form as

Qa q=1 θ q a (µ)A q N u N (µ) =   Q f q=1 θ q f (µ)F q N   , (10)
where (µ N (µ)) j = u Nj , (A q N ) ij = a q ξ pr j , ξ pr i and (F q N ) i = f q (ξ pr i ) .

(11) The output can be expressed as

s N (µ) = Q q=1 θ q (µ) q (u N (µ)) , (12)
or in a vector form

s N (µ) = Q q=1 θ q (µ)L q T N u N (µ) , (13)
where (L q N ) i = q (ξ pr i ). The offline/online decomposition is clear. During the offline step we compute basis functions u N (µ), then form matrices A q N , vectors F q N and vectors L q N . During the online step, for a given parameter µ, we assemble the matrix A N (µ) = Qa q=1 θ q a (µ)A q N , and the vectors

F N (µ) = Q f q=1 θ q f (µ)F q N and L N (µ) = Q q=1 θ q (µ)L q N . We solve the system A N (µ) u N (µ) = F N (µ) , (14
) and finally we can evaluate the output as

s N (µ) = L T N (µ) u N (µ) . (15)
Up to now, we develop a primal-only approach to evaluate the output of interest, but we have not the quadratic convergence effect for outputs, except for compliant cases -a is symmetric and = f . -We introduce now the primal-dual approach for non-compliant cases. The dual problem associated with the output of interest consists in finding Ψ(µ) ∈ X(Ω) such that

a (v, Ψ(µ); µ) = -(v; µ) . (16)
Ψ is denoted adjoint or dual field -note that in the compliant case Ψ = -u. -As in the primal case, we introduce the set S Ψ N that contains evaluation of adjoint for each µ in the set S N

S Ψ N = {Ψ N (µ i ) , ∀µ i ∈ S N } . (17)
We apply the Gram-Schmidt process with respect to the (., .) X inner product to elements of the set S Ψ N , the result are mutually (., .) X -orthonormal basis functions ξ du n , 1 ≤ n ≤ N . The reduced basis space W N du is then defined as

W N du = span{ξ du n , 1 ≤ n ≤ N } .
(18) Note that we have chosen the segregated approach -we could have chosen the integrated one where we have a unique reduced basis space holding both primal and dual basis functions. -In that case, the evaluation of the output is given by s

N (µ) = (u N (µ)) -r pr (ψ N ; µ) . (19) with r pr (v; µ) = f (v; µ) -a(u N (µ), v; µ) and r du (v; µ) = -(v) -a(v, ψ N (µ); µ) .
(20) So now to have to reduced basis approximation, we need to solve the primal and a dual prolem :

a(u N (µ), v; µ) = f (v; µ) ∀v ∈ W Npr and a(v, ψ N (µ); µ) = -(v; µ) ∀v ∈ W N du . (21)
1.2. A posteriori error estimation for elliptic linear problems with affine dependence in parameters

Thanks to (19) we can rapidly compute an estimation for the output of interest s(µ). In this section we introduce an a posteriori error bound that allow us to know if this output estimate is a good enough approximation of the output of interest (see for example [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF], Nguyen et al., 2009]).

Ingredients

We first introduce a positive lower bound α N LB (µ) for α N (µ) for all µ in D, where α N (µ) is the finite elements coercivity constant defined as

α N (µ) = inf w∈N a(w, w; µ) ||w|| 2 X ∀µ ∈ D. ( 22
)
So we can write

0 ≤ α N LB (µ) ≤ α N (µ) ∀µ ∈ D, (23)
where the online computational time to evaluate µ → α N lb (µ) has to be independent of N to compute efficiency the error bound that we introduce here. The successive constraints method (see [START_REF] Huynh | A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants[END_REF]) determine efficiency the lower bound α N LB (µ). Note that in this paper we introduce the SCM because it is implemented in Feel++. There exists other ways to determine α N LB (µ), for example by inspection.

Successive Constraints Method

First we recall that from (3) we assume that the parametrized bilinear form a(., .; µ) depend affinely in µ and can be expressed as

a(u(µ), v; µ) = Qa q=1 θ q a (µ)a q (u(µ), v) ∀u(µ), v ∈ X(Ω), ∀µ ∈ D . (24)
In order to have an offline/online strategy, we will reformulate the expression of the coercivity constant ( 22) by introducing the objective function J obj : D × R Qa → R given by J obj (µ; y) = Qa q=1 θ q (µ) y q where y ∈ R Qa .

Then the FE coercivity constant can be defined by

α N (µ) = inf y∈Y J obj (µ; y) , (25)
where the set Y ∈ R Qa is defined by

Y = y ∈ R Qa | ∃w y ∈ X N s.t. y q = a q (w y , w y ) ||w y || 2 X , 1 ≤ q ≤ Q a . ( 26
)
We also introduce C J as the "coercivity constraint" sample

C J = µ 1 ∈ D, ..., µ J ∈ D . ( 27
)
Lower bound. To determine the lower bound α N LB (µ) we need to define the "continuity constraint" box

B = Qa q=1 inf w∈X N a q (w, w) ||w|| 2 X , sup w∈X N a q (w, w) ||w|| 2 X . ( 28
)
Let C M,µ J the set of M (≥ 1) points in C J closest in the euclidean norm to a given µ ∈ D. The lower bound of the coercivity constant α N (µ) is defined as

α N LB (µ; C J , M ) = min y∈Y LB (µ,C J ,M ) J obj (µ; y) , (29)
where the "lower bound" set Y LB (µ; C J , M ) ∈ R Qa is defined as

Y LB (µ; C J , M ) ≡ y ∈ R Qa | y ∈ B, Qa q=1 β q (µ )y q ≥ α N (µ ) , ∀µ ∈ C M,µ J . ( 30
)
The lower bound defined in ( 29) is a linear optimisation problem, or Linear Program (LP). We observe that the LP (29) contains Q design variables and 2 Q + M inequality constraints. It is important to note that for a given B and the set {α N (µ ) | µ ∈ C J }, the operation count to evaluate µ → α N LB (µ) is independent of N . In actual practice, we have developed a more effcient SCM strategy to build our lower bound, see [START_REF] Vallaghé | A successive constraint method with minimal offline constraints for lower bounds of parametric coercivity constant[END_REF]. Upper bound. We have now a lower bound for the coercivity constraint, but to build the sample C J we also need an upper bound α N U B of this constant. The lower bound of the coercivity constant α N (µ) is defined as

α N U B (µ; C J , M ) = min y∈Y LB (µ,C J ,M ) J obj (µ; y) . (31)
where for given C J , M ∈ N and any µ ∈ D we introduce an " upper bound " set

Y U B (µ; C J , M ) ∈ R Qa as Y U B (µ; C J , M ) ≡ arg inf y∈Y J obj (µ; y) | µ ∈ C M,µ J . ( 32
)
We note that to evaluate µ → α N U B (µ) is independent of N .

Construction of the set C J . Thanks to the previous ingredients, we construct now the set C J using an offline greedy algorithm. First we require a sample Ξ train = {µ 1 , ..., µ ntrain } ⊂ D of n train parameters. A tolerance ∈ (0, 1) is also required to control the error in the lower bound prediction. Start by taking J = 1 and choosing C 1 = µ 1 arbitrarily. Then from J = 2 , the J th parameter µ J selected maximizes the gap between the lower bound and the upper bound of the coercivity constant. That is to say that we perform the algorithm 1.

Algorithm 1 Offline greedy algorithm

while max µ∈Ξtrain α N U B (µ; C J , M ) -α N LB (µ; C J , M ) α N U B (µ; C J , M ) > do µ J+1 = arg max µ∈Ξtrain α N U B (µ; C J , M ) -α N LB (µ; C J , M ) α N U B (µ; C J , M ) C J+1 ← C J ∪ µ J+1 J ← J + 1 end while
At each iteration we add the parameter µ ∈ D that have the worse lower bound approximation to the coercivity constraint sample. As for each µ ∈ C J we have α N U B (µ; C J , M ) = α N LB (µ; C J , M ) it follows from continuity considerations that, for sufficiently large number of iterations, the error tolerance is reached.

Offline/Online strategy of the SCM. This method has obviously an offline/online strategy. During the offline step are computed eigenvalues which serve as bounds of the "continuity constraint" box (28). To determine the upper bound α N U B (µ), the vector y ∈ R Qa , element of the set Y (see 26), is built during the offline step. Then for a given µ ∈ D we determine α N U B (µ) via (31). To apply constraints needed fo the construction of the set Y LB (see 30) we compute eigenvalues associated to parameters µ during the offline step.

A posteriori error estimators

Now that we have a lower bound for the coercivity constant, we introduce the dual norm of the primal (resp. dual) residual Npr (µ) ( resp. N du (µ) ), defined as

Npr (µ) ≡ sup v∈X N r pr (v; µ) ||v|| X = ||ê pr (µ)|| X and N du (µ) ≡ sup v∈X N r du (v; µ) ||v|| X = ||ê du (µ)|| X . (33)
In addition to the dual norm of residuals, equation ( 33) also introduced the Riesz representation of primal (resp. dual) residual : êpr (µ) ( resp. êdu (µ)). Now we define the a posteriori error estimation in terms of the dual norm of residuals and the lower bound for the coercivity constant. In particular for all µ ∈ D and all N we have

|s N (µ) -s N (µ)| ≤ ∆ s N (µ) , (34)
where the a posteriori error estimator on the output ∆ N (µ) is given by

∆ s N (µ) = Npr (µ) 2 α N lb (µ) N du (µ) 2 α N lb (µ) . (35)
Note that (35) can be seen as the product of error estimators on primal and dual solutions.

Offline / online strategy

We now turn to the description of the offline/online strategy for the a posteriori error estimator introduced by (35). We start with the dual norm of the residuals, starting by recalling the expression of the primal residual

r pr (v; µ) = Q f q=1 θ q f (µ)f q (v) - Qa q=1 N j=1 θ q a (µ)u N j (µ)a q (ξ pr j , v), ∀v ∈ X N . ( 36
)
The Riesz representation êpr (µ) verifies

(ê pr (µ), v) X = Q f q=1 θ q f (µ)f q (v) - Qa q=1 N j=1 θ q a (µ)u N j (µ)a q (ξ pr j , v), ∀v ∈ X N (37)
and thus, using linear superposition, it reads

êpr (µ) = Q f q=1 θ q f (µ)Γ q Npr + Qa q=1 N j=1 θ q a (µ)u N j (µ)Υ qj Npr , (38)
where

(Γ q Npr , v) X = f q (v) ∀v ∈ X N , 1 ≤ q ≤ Q f , (Υ qj Npr , v) X = -a q (ξ pr j , v) ∀v ∈ X N , 1 ≤ q ≤ Q a , 1 ≤ j ≤ N . (39)
Consequently we have

Npr (µ) 2 = C f f Npr (µ) + 2 N j=1 u N j C f a Nprj (µ) + N j=1 N j =1 u N j u N j C aa Nprjj (µ) , (40)
where for 1 ≤ j, j ≤ N :

C f f Npr (µ) = Q f q=1 Q f q =1 θ q f (µ) θ q f (µ) Φ qq Npr , C f a Nprj (µ) = Qa q=1 Q f q =1 θ q a (µ) θ q f (µ) Φ qjq Npr , C aa Nprjj (µ) = Qa q=1 Qa q =1 θ a q (µ) θ a q (µ) Φ qjq j Npr . (41)
Here

Φ qq Npr ( 1 ≤ q, q ≤ Q f ) , Φ qjq Npr ( 1 ≤ q ≤ Q a ,1 ≤ j ≤ N , 1 ≤ q ≤ Q f ) and Φ qjq j

Npr

( 1 ≤ q, q ≤ Q a ,1 ≤ j, j ≤ N ) do not depend on parameter µ and are defined as

Φ qq Npr = (Γ q Npr , Γ q Npr ) X , Φ qjq Npr = (Υ qj Npr , Γ q Npr ) X and Φ qjq j Npr = (Υ qj Npr , Υ q j Npr ) X . (42)
Now let us consider the expression of the dual residual

r du (v; µ) = - Q l q=1 θ q l (µ)l q (v) - Qa q=1 N j=1 θ q a (µ)Ψ N j (µ)a q (v, ξ du j ), ∀v ∈ X N . ( 43
)
The Riesz representation êdu (µ) reads

êdu (µ) = Q l q=1 θ q l (µ)Γ q N du + Qa q=1 N j=1 θ q a (µ)Ψ N j (µ)Υ qj N du , (44)
where

(Γ q N du , v) X = -l q (v) ∀v ∈ X N , 1 ≤ q ≤ Q l , (Υ qj N du , v) X = -a q (v, ξ du j ) ∀v ∈ X N , 1 ≤ q ≤ Q a , 1 ≤ j ≤ N . (45)
Consequently êdu (µ) 2 can be written as

N du (µ) 2 = C ll N du (µ) + 2 N j=1 Ψ N j C la N du j (µ) + N j=1 N j =1 Ψ N j Ψ N j C aa N du jj (µ) , (46)
where for 1 ≤ j, j ≤ N :

C f f N du (µ) = Q f q=1 Q f q =1 θ q f (µ) θ q f (µ) Φ qq N du , C f a N du j (µ) = Qa q=1 Q f q =1 θ q a (µ) θ q f (µ) Φ qjq N du , C aa N du jj (µ) = Qa q=1 Qa q =1
θ a q (µ) θ a q (µ) Φ qjq j N du .

(47)

Here Φ qq N du ( 1 ≤ q, q ≤ Q f ) , Φ qjq N du ( 1 ≤ q ≤ Q a ,1 ≤ j ≤ N , 1 ≤ q ≤ Q f ) and Φ qjq j N du
( 1 ≤ q, q ≤ Q a ,1 ≤ j, j ≤ N ) do not depend on parameter µ and are defined as

Φ qq N du = (Γ q N du , Γ q N du ) X , Φ qjq N du = (Υ qj N du , Γ q N du ) X and Φ qjq j N du = (Υ qj N du , Υ q j N du ) X . (48)
During the offline stage, we compute all µ-independent quantities and we store them in a database. During the online step, for a given µ ∈ D, we evaluate the µ-dependent terms and assemble the residual norm terms using the online solution.

Empirical Interpolation Method

Finally we describe the Empirical Interpolation Method (see [START_REF] Barrault | An empirical interpolation method: application to efficient reduced-basis discretization of partial differential equations[END_REF], [START_REF] Grepl | Efficient reduced-basis treatment of nonaffine and nonlinear partial differential equations[END_REF]). We consider the non-linear µ-dependent function of sufficient regularity g (µ; x; u(µ)) which is non-affine in parameters and depends on the solution of a parametrized PDE u(µ). We are interested in approximating g (µ; x; u(µ)) by a reduced-basis expansion g M (µ; x; u(µ)) which is affine in the parameters such that

g M (µ; x; u(µ)) = M m=1 β m (µ; u(µ))q m (x) . ( 49
)
To this end we introduce a nested sample set

S M = {µ 1 , ..., µ M } ∈ D M and associated nested space W M = span{ξ m ≡ g (µ m ; x; u(µ m )) , 1 ≤ m ≤ M } in which the approximation g M (µ; x; u(µ)) shall reside.
We first introduce Ξ a suitably large but finite-dimensional parameter set in D. The first sample point µ 1 is picked into Ξ and assuming that ξ 1 = 0, we define :

S M = {µ 1 }, ξ 1 ≡ g (µ 1 ; x; u(µ 1 )) and W M = span{ξ 1 } .
For M 2, we determine µ M from a Greedy algorithm and deduce the associated ξ M :

µ M = arg max µ∈Ξ inf z∈W M -1 ||g(µ; .; .) -z|| L ∞ (Ω) , and ξ M = g (µ M ; x; u(µ M )) (50) from which we complete S M = S M -1 ∪ {µ M } and W M = W M -1 ⊕ span{ξ M }.
The coefficients β m of the particular linear combination (49) are determined through interpolation points

t 1 , . . . , t M ∈ Ω such that M m=1 β m (µ; u(µ))q m (t i ) = g(µ; t i ; u(µ)) ∀ t i (i = 1 . . . M ) . ( 51
)
The first interpolation point t 1 is chosen such that first basis function ξ 1 is maximum, and q 1 consists in the normalization of ξ 1 :

t 1 = arg sup x∈Ω |ξ 1 (x)|, q 1 = ξ 1 (x) ξ 1 (t 1 ) and B 1 11 = q 1 (t 1 ) = 1 . ( 52
)
For M 2, we look for the vector σ

M -1 = σ M -1 i i=1,...,M -1 obtained through interpolation points M -1 j=1 σ M -1 j q j (t i ) = ξ M (t i ) , 1 ≤ i ≤ M -1 , (53)
and then we compute the residual

r M (x) = ξ M (x) - M -1 j=1 σ M -1 j q j (x) . ( 54
)
which allows to compute the next interpolation point t M together with the basis functions q M (x) for M 2 along with the interpolation matrix

B M t M = arg sup x∈Ω |r M (x)|, q M (x) = r M (x) r M (t M ) and B M ij = q j (t i ) , 1 ≤ i, j ≤ M . ( 55
)
Once all the interpolation points t M and the basis functions q M have been computed (offline part), the computation of the approximation g M (µ; x; u(µ)) for a given µ (online part) consists in finding the coefficients β m (µ, u(µ)) of the linear combination by solving :

M j=1 B M ij β j (µ; u(µ)) = g (µ; t i ; u(µ)) , 1 ≤ i ≤ M . (56)
In other words, coefficients of EIM expansion are solutions of

     q 1 (t 1 ) 0 • • • 0 . . . . . . . . . . . . q 1 (t M -1 ) • • • q M -1 (t M1 ) 0 q 1 (t M ) • • • • • • q M (t M )           β 1 . . . β M -1 β M      =      g (µ; t 1 ; u(µ)) . . . g (µ; t M -1 ; u(µ)) g (µ; t M ; u(µ))      . ( 57
)
Remark 1. In the most general case, the function g depends not only on x but also the solution u(µ) of the problem. During the online step, we need to evaluate quickly the function g at interpolation points (see 57), it implies to have a rapid evaluation of u(µ) which is in fact replaced by u N (µ) and we have precomputed the reduced basis functions associated to u N (µ) at the interpolation points (t i ) i=1,...,M .

1.4. Elliptic non-affinely parametrized non-linear equations

1.4.1. General problem settings
We consider the following problem : given µ ∈ D ⊂ R p , evaluate the output of interest :

s(µ) = (u(µ)) , (58)
where u(µ) ∈ X(Ω) satisfies g(u(µ), v; µ) = 0 ∀v ∈ X(Ω) .

(59) Here we consider that ( 59) is a non-linear system of N equs equations. To deal with (59), the Newton algorithm is used to find zero of g(u(µ), v; µ). Let k u(µ) the solution at the k th iteration of Newton algorithm, g k u(µ), v; µ the function g applied to the solution k u(µ) and j u(µ), v; µ; k u(µ) the jacobian of g(u(µ), v; µ) applied to the solution k u(µ). It reads

j u(µ), v; µ; k u(µ) =         ∂g 1 u(µ), v; µ; k u(µ) ∂u 1 • • • ∂g 1 u(µ), v; µ; k u(µ) ∂u Nequs . . . ∂g Nequs u(µ), v; µ; k u(µ) ∂u 1 • • • ∂g Nequs u(µ), v; µ; k u(µ) ∂u Nequs         . ( 60
)
Now, we make the assumption for the offline/online procedure by assuming that we can approximate g (., .; µ), j ., .; µ; k u(µ) and the linear form (.; µ) respectively by g AD (., .; µ) , j AD ., .; µ; k u(µ) and AD (.; µ) which are functions that are affine in parameters µ. That is to say that for given Q g , Q j and Q the EIM determines (M g q ) q=1,..,Qg , (M j q ) q=1,..,Qj and (M q ) q=1,..,Q such that we can write :

g AD k u(µ), v; µ = Qg q=1 M g q m=1 β qm g (µ; k u(µ)) g qm (v) , (61)
j AD u(µ), v; µ; k u(µ) = Qj q=1 M j q m=1 β qm j (µ; k u(µ)) j qm (u(µ), v) , (62)
and

AD (v; µ) = Q q=1 M q m=1 β qm (µ) qm (v) . ( 63
)
The Newton method is an iterative method (see algorithm 2). Starting with an initial guess which is reasonably close to the true root the method consists, for each step, in solving the system of linear equations

j AD u(µ), v; µ; k u(µ) δ k u(µ) = -g AD k u(µ), v; µ , (64)
for the increment δ k u(µ) defined by

δ k u(µ) = k+1 u(µ) -k u(µ) . (65)
Algorithm 2 u(µ) = N ewton(g( k u(µ), v; µ), j(u(µ), v; µ; k u(µ)), initial guess, tol)

k ← 1 1 u(µ) ← initial guess e ← ∞ while e > tol do solve j u(µ), v; µ; k u(µ) δ k u(µ) = -g k u(µ), v; µ e ← k+1 u(µ)-k u(µ) k u(µ) k ← k + 1 end while u(µ) ← k+1 u(µ) 1.4.2

. Reduced basis method

The finite element discretization of problem ( 58)-( 59) is as follows : given a µ ∈ D, evaluate

s N (u N (µ)) = AD (u N (µ); µ) , (66)
where

u N (µ) ∈ X N satisfies g AD (u N (µ), v; µ) = 0 ∀v ∈ X N .
(67) Here is mentioned g AD (., .; µ) instead of g(., .; µ) to recall that we use approximations described in (61),( 62) and ( 63). The reduced basis solution u N (µ) can be expressed as

u N (µ) = N i=1 u N i (µ) ξ i . ( 68
)
Considering the k th iteration of the Newton algorithm, by taking v = ξ i for i = 1, ..., N and using (68) we can write

N n=1 j AD ξ n , ξ i ; µ; k u N (µ) δ k u N n (µ) = -g AD k u N (µ), ξ i ; µ , (69)
where δ k u N n (µ) is the n th component of the increment, defined as :

δ k u N n (µ) = k+1 u N n -k u N n
. Now we use the affine parameter dependence to construct an efficient offline/online strategy. Choosing test functions as v = ξ i , i = 1, . . . N , the equation ( 69) can be expressed as

N n=1   Qj q=1 M j q m=1 β qm j (µ; k u N (µ)) j qm (ξ n , ξ i )   δ k u N n (µ) = - Qg q=1 M g q m=1 β qm g µ; k u N (µ) g qm (ξ i ) . ( 70
)
The equivalent matrix form is

  Qj q=1 M j q m=1 β qm j (µ; k u N (µ)) J qm N   k+1 u N (µ) -k u N (µ) = - Qg q=1 M g q m=1 β qm g µ; k u N (µ) G qm N . (71)
The

unknown k u N (µ) ∈ R N is defined as k u N (µ) n = k u N n (µ). The matrix J qm N ∈ R N ×N and the vector G qm N ∈ R N are defined as (J qm N ) i,n = j qm (ξ n , ξ i ) and (G qm N ) i = g qm (ξ i ), 1 i, n N. (72)
The output can be expressed as

s N (µ) =   Q q=1 M q m=1 β qm (µ) qm (u N (µ))   , (73)
or in a vector form

s N (µ) =   Q q=1 M q m=1 β qm (µ) L qm T N   u N (µ) , (74)
where (L qm N ) i = qm (ξ i ). During the offline step we compute basis functions u N (µ), then form matrices J qm N , vectors G qm N and vectors L qm N . During the online step, for a given parameter µ and a given solution k u N (µ) we update the jacobian matrix

J N µ; k u N (µ) = Qj q=1 M j q m=1 β qm j (µ; k u N (µ))J qm N , and the vector G N (µ) = Qg q=1 M g q m=1 β qm g (µ; k u N (µ))G qm N , to find k+1 u N (µ) such that J N µ; k u N (µ) k+1 u N (µ) -k u N (µ) = -G N µ; k u N (µ) , (75)
until convergence. Then, when we have the solution u N (µ) we assemble the vector

L N (µ) = Q q=1 M q m=1 β qm (µ)L qm N , (76)
and finally we can evaluate the output as

s N (µ) = L T N (µ) u N (µ) . ( 77
)

Computational framework for reduced basis

We now describe the RB framework used is the Feel++ RB framework. Its relies on Feel++ which we describe briefly before turning to the RB layer.

Feel++ : principles and design

The library Feel++ 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 and in our case the DSEL splits low level mathematics and computer science on one side leaving the Feel++ developer to enhance them and high level mathematics as well as physical applications to the other side which are left to the Feel++ user. This 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, yet with a simple and seamless interface, tools such as interpolation or the clear translation of a wide range of variational formulations into the variational embedded language. Combined with this robust engine, lie also state of the art arbitrary order finite elements -including handling high order geometrical approximations, -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 the numerical approximations.

The code 1 illustrates the clear and easy implementation -building of mesh and function spaces, writing of variational formulation -of a laplacian problem with homogeneous Dirichlet conditions, provided by the Feel++ library : // post processing auto e = exporter ( _mesh = mesh ); e -> add ( " u " , u ); e -> save (); } As to build the reduced basis during the offline step of the RBM we need to use the FE discretization, we used a recent development of Feel++ which allows to operate on parallel computers. In order to create a parallel computing code, we use some strategies of domain decomposition methods with the MPI technology. A feature of our library is that all MPI communications are seamless, thanks to DSEL. Thus, the library Feel++ provides a parallel data framework : we start with automatic mesh partitioning using Gmsh [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF](Chaco/Metis) -adding information about ghosts cells with communication between neighbor partitions; -then Feel++ data structures are parallel such as meshes, (elements of) function spaces -parallel degrees of freedom table with local and global views; -and finally Feel++ uses the library PETSc [START_REF] Balay | PETSc Web page[END_REF], Balay et al., 2012a, Balay et al., 1997] which provides access to a Krylov subspace solvers(KSP) coupled with PETSc preconditioners such as Block-Jacobi, ASM, GASM. The last preconditioner is an additive variant of the Schwarz alternating method for the case of many subregion, see [START_REF] Smith | Domain decomposition: parallel multilevel methods for elliptic partial differential equations[END_REF]. A complete description of this HPC part with some blood flow applications is done in the thesis [Chabannes, 2013]. The RB framework, depicted in figure 1, provides an interface to the reduced basis methodologies presented so far and automatically generated for the User different instantiations of the reduced basis applications namely command line executable as well as Python and Octave interfaces. These interfaces design follow the simple input-output relationships:

Feel++ reduced basis framework

[s 1 (µ), s 2 (µ), . . . , s O (µ)] = (µ) (78)

where is now the simulation software which takes the input parameter set µ and s i are the outputs of the software which is a set of performance metrics and, if available, the associated error bounds. The User needs to provide the specifications of the model: parameter space, geometry, variational formation with affine decomposition (possibly using EIM).

The design of the C++ classes RB framework is illustrated by the figure 2. The class ParameterSpace generates and stores parameter samplings which are then used by CRBSCM and EIM. Various sampling strategies are available. CRBModel is the interface that the User model must support to use the RB framework and this is through this interface that the User provides the model specifications described previously.

The class CRB represents the certified reduced basis implementation for elliptic and parabolic. For elliptic problems by default the Greedy algorithm is used. For parabolic problems, the POD/Greedy algorithm is used ; POD in time and Greedy in space. Note that in both cases, the user can specified in a file the sampling S N (used to build the space W Npr and W N du ). S N can also contains log-equidistributed or log-randomized parameters. The class CRBDB handles the database storage using the Boost.serialization library. Concerning the online step, dense matrices and vectors are manipulated via the library Eigen [START_REF] Guennebaud | Eigen v[END_REF]. For linear equations linear solvers from Eigen are used. Non-linear solvers from PETSc are used to deal with non-linear equations. In order to interact with PETSc solvers, Eigen::Map<> is the bridge to communicate the non-linear data between Eigen and PETSc. The class CRBSCM implements the SCM algorithm thanks to standard and generalized Eigensolvers from SLEPc [START_REF] Hernandez | SLEPc: A scalable and flexible toolkit for the solution of eigenvalue problems[END_REF]. The class EIM implements the EIM algorithm. To build the reduced basis, the offline step of the method can be very expensive. Once this expensive offline step is done, we save all scalar product results from the projection of matrices or vectors on the reduced basis in a database. Hence we can reuse an existing basis to perform online computations. We store the projection of matrices and vectors on the reduced basis, see 11. To estimate the error, we store scalar products defined in ( 42) and (48). For the SCM, all quantities computed during the offline step mentioned in section 1.2.2 are stored, they refer to equations ( 28), ( 26) and (30). Concerning EIM, the interpolation matrix and basis functions are stored, see ( 52) and (55). To save objects in a database, we use the serialization process provided by the Boost.serialization library. The framework is designed so that enriching an existing basis is possible.

A sensitivity analysis can be performed using the scientific library OpenTURNS [START_REF] Works | OpenTURNS Web page[END_REF], specialized in the treatment of uncertainties. Through python scripts, the User set a range and a distribution for the inputs from which OpenTURNS build a sampling, given to CRBApp. Running the online step of the RBM, CRBApp supplies the associated set of outputs to OpenTURNS which can then compute quantities of interest such as standard deviation, quantiles and Sobol indices.

If the User wants to test the RBM on a model, CRBApp will generates a log-randomized or log-equidistributed sampling of parameters on which the output will be evaluated. To verify the RB approximation of the output, CRBApp provides also the output evaluated by using the finite element discretization. By assuming than, for a given parameter µ ∈ D, the solution field u N (µ) obtained using the finite element discretization is the true solution, the error on the solution filed is then computed in L 2 and H 1 norm. In other words the class provides

e L2 = ||u N (µ)-u N (µ)|| L 2 (Ω) ||u N (µ)|| L 2 (Ω)
and e H1 =

||u N (µ)-u N (µ)|| H 1 (Ω) ||u N (µ)|| H 1 (Ω)
. It also computes the error on the output :

e output = |s N (µ)-s N (µ)| |s N (µ)|
.

Parallel strategy

The RB framework supports parallel architectures by relying on the Feel++ parallel data structures. All data associated to the reduced basis (scalars, dense vectors and matrices, parameter space samplings) are actually 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. Finally, in the implementation of the EIM algorithm, we need to compute the arg max over a fine sampling of infinity norm of a quantity that lives on the mesh, see equation ( 50). Each processor computes its local contribution and the a "max" MPI reduction is done.

Practical use of EIM

We consider now a practical example of the usage EIM taken from the non-linear electro-thermal problem described in section 3.1. The associated weak formulation can be expressed as :

find U = (V, T ) ∈ X h s.t. ∀V ∈ X h and ∀µ ∈ D : g(U , V ; µ) = 0 . ( 79
)
where V and T are respectively the electrical potential and the temperature, and X h is a composite space (product of spaces X V h and X T h associated with V and T ). In this model, g depend on a term σ(T ), which is non-affine in parameter µ :

σ(T ) = µ 0 1 + µ 1 (T -T 0 ) , ( 80
)
where T 0 is a set constant. As a consequence, g is non-affine in parameter µ and the associated affine decomposition can not be obtained without an EIM decomposition of the term σ(T ). To recover an affine decomposition of ( 80), the User has to use EIM class by instantiating an EIM object as shown in code 2. This EIM object gives a direct access to the basis functions q M and the coefficients β M of the EIM expansion (described in 1.3). Now, let us see how to have an expansion of a function which does not depend on the solution of the model. Say we want the expansion of g(µ; x) = sin (µ 0 π x) where µ 0 is a parameter of the model. Here there is no need to compute the solution of the model to have the expansion of g(µ; x). So we use the keyword eim_no_solve, like illustrated in code 3.

Listing 2. EIM approximation of a non-affine parameter and solution dependent term

Listing 3. EIM approximation with no dependence on solution // eim expansion of sin (µ0 π x) auto eim_sin = eim ( _model = eim_no_solve ( FEMModel ) ) , _element =T , _parameter = mu , // µ _expr = sin ( ref ( mu (0) )* pi * Px () ) , _space =X h , _name = " eim_sin " , _sampling = Pset );

Numerical experiments

The Laboratoire National des Champs Magnétiques Intenses (LNCMI) is a French large scale facility ( [START_REF] Debray | The grenoble high magnetic field laboratory as a user facility[END_REF]) enabling researchers to perform experiments in the highest possible magnetic field (up to 35 T static field). High magnetic fields are obtained by using water cooled resistive magnets (cf. fig. 3 and4) connected with a 24 MW power supply (see [F. Debray et al., 2012] and references included for details).

The design and optimization of these magnets require from an engineering point of view the prediction of "quantities of interest," or performance metrics, which we shall denote outputs -namely magnetic field in the center, maximum stresses, maximum and average temperatures. These outputs are expressed as functionals of field variables associated with a set of coupled parametrized PDEs which describe the physical behavior of our magnets. The parameters, which we shall denote inputs, represent characterization variables -such as physical properties, heat transfer coefficients, water temperature and flowrate, and geometric variables in optimisation studies. To evaluate these implicit input-output relationships, solutions of a multi-physics model involving electro-thermal, magnetostatics, electro-thermal-mechanical and thermo-hydraulics are requested. It should be noted that this model is non-linear as the material properties depend on temperature. In practice these evaluations represent a huge computational time but they are mandatory to improve the magnet design as we can no longer rely on common physical sense.

To significantly reduce this computational time, we choose to use model order reduction strategies, and specifically to use the reduced basis method presented in the previous sections which is well adapted to tackle this question. First, we focused on the electro-thermal behaviour of the resistive magnets. We present now the non-linear electro-thermal coupled problem developed for this purpose. Design issues and challenges attached to this application are then explicited. Finally the Feel++ reduced basis framework is applied on two examples directly connected with actual LNCMI developments.

A non-linear electro-thermal model for High Fields Magnets

Two main technologies are developed at the LNCMI for the resistive magnets, namely the Bitter and Polyhelix techniques. A typical 24 MW magnet will consist in a set of Polyhelices and Bitter magnets (or inserts in LNCMI standard terminology) powered separatly by a 12 MW power supply. Polyhelices insert are in turns made of concentric copper alloy tubes, in which a helix have been cut by spark erosion techniques, electricaly connected in series. The helix cut in each tube may be filled with some glue or let free. In this case some insulator are introduced periodicaly to prevent electrical contact between turns.

Figure 3. A Polyhelix insert. The helices are actually copper alloy tubes cut by a spark erosion technique. Epoxy glue may be introduced into the slit to ensure the electrical insulation between turns : this is the so-called "Longitudinally cooled helices" as the water flow is "longitudinal". An another possibility to ensure electrical insulation is to introduce periodically some insulators into the slit : this is the so-called "Radially cooled helices" as the water may flow from the inner to the outer radius througth the open slit. Small holes are drilled to enable the water flow throught the magnet. Some tie-rods are inserted to ensure a good electrical contact between each plate by applying a pre-compression.

The applied electrical current is about 30 kA in each insert. This leads to important Joules losses within the insert. A water flow (about 150 l/s) is then necessary to cool down the insert and to avoid that the temperature reaches some given threshold.

The temperature T in the magnets is given by the following electro-thermal coupled model, where V stands for the electrical potential :

-∇ • (σ(T ) ∇V ) = 0, -∇ • (k(T ) ∇T ) = σ(T ) ∇V • ∇V. (81)
σ(T ) and k(T ) (respectively the electrical and thermal conductivity of the material) depend on T , hence the non linearity of this model :

σ(T ) = σ 0 1 + α(T -T 0 ) and k(T ) = σ(T )LT (82)
where σ 0 denotes the electrical conductivity at T = T 0 . α and L are characteristics of material, namely the ratio resistivity-temperature and the Lorentz number.

Dirichlet boundary conditions on V are applied to the surfaces connecting the magnet to the power bus bar. On all other boundary Neumann conditions are applied, as no current flows out the magnets apart from the top and bottom regions : V = 0 on top and V = V D on the "connection" surfaces, -σ(T ) ∇V • n = 0 on other boundaries.

The applied V D will be estimated such that the total electrical current in the magnet is equal to some given value.

We can assume that there is no thermal exchange on the "connection" surfaces, nor on the interfaces with insulators. The thermal exchange between copper and cooling water is modelled by a forced convection conditions (with h the heat transfer coefficient, and T w the water temperature) :

-k(T ) ∇T • n = 0 on the "connection" surfaces, -k(T ) ∇T • n = h(T -T w ) on regions in contact with water.

h is obtained by mean of classical thermohydraulics correlations given the water flowrate.

Design issues and Challenges

In the international race for higher magnetic field we need to push our magnet technologies to their limits. The actual design of our inserts are mainly limited by thermal constraints for the inner parts and by mechanical constraints for the outer parts. Two routes are actively investigated to go beyond these limits : one consists in developing the "Radially cooled helices" (see comment on fig 3) which are less sensitive to thermal limitations from their design principles; the other consists in looking for materials with improved mechanical properties. This last aspect is a more long term research.

LNCMI current projects involves the use of those "Radially cooled helices" to achieve higher magnetic field by breaking the thermal limits. From a designer point of view it requires to control the temperature, more precisely the average and standard deviation temperature, in the magnets for every operating conditions (ie. applied potential V D and cooling conditions T w , h). Furthermore the material properties are not precisely known : the copper producer only gives upper and lower bounds for σ 0 ; as for α and L they are only known from litterature to be in a given interval. Accounting for these uncertainties is also important to be confidant in our design. This is quite challenging as the temperature field is solution of the non-linear system of PDEs (81).

To reach this goal we need fast, accurate and reliable estimates for the temperature field. This need may be adressed by controlling a posteriori error for T associated with an anisotropic mesh adaptation strategy ( [START_REF] Prud'homme | Advances in feel++ : A domain specific embedded language in C++ for partial differential equations[END_REF]). This approach is however not sufficient nor efficient in terms of required computational ressources. An efficent alternative is to use the Reduced basis method (see 1.4.2) to greatly reduce the computional cost. In our case this method has to be combined with the Empirical Interpolation Method (see 1.3) to deal with non-affine dependence parameters, arising from the righthand side of ( 81) and ( 82). This new approach implemented using Feel++ reduced basis framework makes it possible to carry out parametric studies and sensitivity analysis at reasonnable cost which are mandatory to improve the design. In the sequel we will consider the input parameter µ ∈ R 6 defined by : µ = (σ 0 , α, L, V D , h, T w ) (83)

A first example of a parametric study on a small geometry is presented in 3.4. It also serves as a validation of the framework. A second example is given in 3.5 for a larger problem, illustrating the parallel computing possibilities of Feel++ RB framework. This second example also shows the computation of quantities of interest for sensitivity analysis -as a first step toward uncertainties quantification-adressed to applications such that design under uncertainties.

Reduced basis approximation of the problem

As the material properties makes our model non-affinely parametrized and non-linear, the use of reduced basis method in such a context involves to proceed as described in section 1.4. The output (equ. ( 58)) we focus on in the next sections is the mean of temperature T (µ) :

s(T (µ)) = 1 |Ω| Ω T (µ)dΩ (84)
The solutions V (µ) and T (µ) of our coupled problem (c.f. section 3.1) are respectively zeros of the two functions g V (V (µ), φ V ; µ) and g T (T (µ), φ T ; µ) expressed as follows :

g V (V (µ), φ V ; µ) = Ω σ(T )∇V • ∇φ V - top∪bottom σ(T )(∇V • n)φ V + top∪bottom σ(T )γ h s V φ V (85) - top∪bottom σ(T )(∇φ V • n)V - top σ(T )γ h s V D φ V + top σ(T )(∇φ V • n)V D g T (T (µ), φ T ; µ) = Ω k(T )∇T • ∇φ T + cooling hT φ T - Ω σ(T )φ T ∇V • ∇V - cooling hT w φ T ( 86
)
Remark 2. To handle seamlessly the Dirichlet conditions in our RB framework we have chosen to treat these conditions weakly using the Nitsche formulation. γ and h s denote respectively a penalisation parameter and h s the characteristic mesh size.

We can now use a Newton algorithm to solve g(< V (µ), T (µ) >, < φ V , φ T >; µ) = g V (V (µ), φ V ; µ) + g T (T (µ), φ T ; µ) = 0 (c.f. section 1.4.1, algorithm 2) where the residual and jacobian read:

r(< V (µ), T (µ) >, < φ V , φ T >; µ) = g V (V (µ), φ V ; µ) + g T (T (µ), φ T ; µ) (87) j(< V (µ), T (µ) >, < φ V , φ T >; µ) =      ∂g V (V (µ), φ V ; µ) ∂V ∂g V (V (µ), φ V ; µ) ∂T ∂g T (T (µ), φ T ; µ) ∂V ∂g T (T (µ), φ T ; µ) ∂T      (88)
The terms σ and k (defined in ( 82)) which appears in g V and g T present a non-affine dependence in parameters. Consequently, we need to use of the Empirical Interpolation Method (c.f. section 1.3) to build the affinely parametrized approximations σ Mσ and k

M k σ Mσ (T (µ), x; µ) = Mσ m=1 β σ m (T (µ); µ)q σ m (x) and k M k (T (µ), x; µ) = M k m=1 β k m (T (µ); µ)q k m (x) (89)
The affine decompositions ( 61)-( 62) can finally be performed in order to apply the offline/online procedure described in 1.4.2.

Parametric study -example for a Bitter magnet

The LNCMI is currently involved in the Hybrid Magnet Project, which consists in the assembly of a resistive inner coil with a superconducting outer one. This technology allows to generate the highest continuous magnetic field for a given electrical power installation (43T for 24MW at LNCMI). The inner coil is a combination of two magnets technology (Bitter and polyhelix), powered independently (see [F. Debray et al., 2012] and [A. Bourquard et al., 2010]).

To reach the target magnetic field, the Bitter magnet has to be redesigned to generate more than 9 Tesla which corresponds to the actual design. This may be achieved by increasing locally the current density. From a pratical point of view it means modifying the stacking of copper plates (see comment on fig 4). This change will also leads to locally increase the temperature.

To determine how much we can possibly increase the current density while keeping a reasonable mean temperature to not induce damage on materials we have carried out a parametric study. This study consists in setting all input parameters except V D . The current density (denoted j c ) and the applied potential V D in the Bitter magnet are related such that

j c = V D r int σ 0 θ (90)
where r int is the inner radius of the disk, and θ corresponds to the angle of the sector.

The input parameter V D can then be chosen such that the current density varies from 30.10 6 to 100.10 6 A.m -2 .

The study was performed with classical FE method by running calculations for a given set of parameters and with CRB method in order to validate the framework. The simulation illustrated above (fig. 5) have been performed on 8 processors (domain divided into 8 partitions), using the solver GMRES and the Additive Schwarz Method (GASM) as preconditioner. In a computational time point of view, the use of reduced basis framework for this computation represents a gain of 2 or 3 orders of magnitude in comparison with the classical FE model.

In the current design the mean temperature is about 40C • . Allowing this value to reach 60C • we can increase the current density by a factor 1.5 as can be seen on fig 6. This allow us to redefine the stacking of the copper plates to safely reach 10 tesla. The new Bitter magnet design will be based on this result.

Towards uncertainty quantification -example for a helix magnet

As stated in introduction to this section, using radially cooled helices for the inner most helices of a polyhelix magnet (see [F. Debray et al., 2012]) allows to reach higher current density and hence to reach higher magnetic fields. This is possible as the radially cooled helices are by concept less thermally limited than longitudinally ones (see comments on fig 3). The water flows from the inner radius to the outer radius of the helix is indeed more efficient to cool it down. However since insulators have to been introduced in between each turn (see fig 7) some hotspots are expected in these regions. Accurate estimates of the temperature are a key element for the design. In particular we want the mean temperature in the helix to remain below some threshold. This allows to limit the maximum temperature reached in the insulators and thus to keep them safe. To do so we must take into account the uncertainties on the copper alloy properties (σ 0 , α and L), and the uncertainties on the operating parameters V D and cooling conditions T w , h. This can be done by performing a sensitivity analysis. The ranges of input parameters indicated in the table below are chosen from literature and experimental results.

Parameter Range σ 0 [10 6 Ω -1 m -1 ] [50; 50, 2] α [10 -3 ] [3, 3; 3, 5] L [10 -8 ] [2, 5; 2, 9] V D [V] [0, 14; 0, 15] h [W K -1 m -2 ] [70000; 90000] T w [K] [293; 313]
Table 1. Input parameters ranges Sensitivity analysis is managed by linking the Feel++ reduced basis framework with the library Openturns, dedicated to the treatment of uncertainties. To perform this study, OpenTurns build a sample of input parameters, determined upon a given probability distribution in the wanted ranges. The RB framework is then used to compute the associated set of outputs, from which the mean and the standard variation are deduced by OpenTurns.

Parallel computing possibilities available within CRB framework allows to apply model order reduction on such a complex geometry (see the left part of fig 8) with 1.5 10 6 dofs. Considering duration of both offline and online steps, we assess that currently the use of reduced basis provides a gain in time compared with FE model over 8 evaluations. As one can see in (89), the coefficients β σ m and β q m which takes part in the EIM decomposition of terms σ and k depend on the solution T (µ). Currently, the computation of EIM online step is still dependent on the dimension of FE discretization space when building the right hand side of the system (57) which requires -currently undertaken -new code design and implementation. We would then expect a gain of 2 or 3 orders of magnitude as for applications presented e.g. in [START_REF] Schenone | High performance computing for the reduced basis method, application to natural convection[END_REF].

Nevertheless, the number of runs required for a sensitivity analysis with OpenTurns is so large (typically 10 000) that we cannot perform it with the computionnal ressources available. For this reason we choose to select a "representative" helix sector, since the temperature behaviour is regular and exhibits some symmetry (see the right part of fig 8).

The following results -on the helix sector -have been obtained from an uniform distribution of the inputs in the ranges defined in table 1, with a sample of size 500. We shall note that -as this example described only the first step towards sensibility analysis-, such a sample is sufficient since we note compute the Sobol indices.

The two following quantities are essential since they give a "reference" value for the mean temperature in the given ranges, associated to the standard deviation which can be considered as a "safety interval" around the previous reference value.

Mean of outputs [K]

368.66 Standard deviation [K] 6.22

Table 2. Mean and standard deviation

The values obtained for mean of outputs and standard deviation ensure that for the given input parameters ranges (table 1), the mean temperature in the magnet will be included in [362, 44; 374, 88] K. The largest obtained value -the most critical-can then be "refined" using the quantiles. Quantiles correspond to the threshold q(γ) we are sure not to exceed with a given probability γ. Considering Y the set of outputs obtained from the sampling, this reads as:

Find q(γ) such that P (Y < q(γ)) > γ (91) For the same configuration as for the mean and standard deviation (table 2), quantiles with respectively γ = 80% and γ = 99% have been computed : Probability (γ) Quantile [K] 80% 371.297 99% 376.014

Table 3. Quantiles These results (table 3) are also of great interest for our magnet control system in order to early detect unexpected thermal behaviours and thus to better predict incidents.

Conclusions and perspectives

Feel++ Reduced Basis framework has been successfully applied to industrial class problems as illustrated in the examples above. It shows its ability to deal with a 3D elliptic non-affinely parametrized non linear problem. From a practical point of view, solving such non-linear coupled problems raise the challenge of management of large databases. The framework currently provides a posteriori error estimators for elliptic linar problems (see section 1.2). Their extension for the certification of the RB approximation in the case of non-affinely parametrized and non-linear problems (see for example [Cuong, 2005]) needs to be implemented. With adequate a posteriori error estimators, a greedy algorithm can then be applied to construct an optimal reduced basis for the electro-thermal coupled problem introduced in section 3.1.

More mid-terms developments involves the implementation of automatic differentiation to greatly help to account for the geometric input parameters. Moreover, the integration of stochastic approach possibilities such as chaos polynomial would simplify the uncertainty quantification. From applications point of view in the context of LNCMI, the electro-thermal model presented here would be enrich with other physics -such as magnetostatic or elasticity -and evolve to handle transient simulations. Finally the possibilities offered by lego simulation linked with domain decomposition techniques (see [START_REF] Vallaghé | The static condensation reduced basis element method for a mixedmean conjugate heat exchanger model[END_REF]) are very appealing for such large applications, involving multi-physics simulations on complex geometries. It would for example allow to perform numerical experiments on complete set of magnets by assembling magnet components.

Listing 1 .

 1 Laplacian problem with homogeneous Dirichlet conditions # include < feel / feel . hpp > int main ( int argc , char ** argv ) { using namespace Feel ; Environment env ( _argc = argc , _argv = argv , _desc = feel_options () , _about = about ( _name = " laplacian " , _author = " Feel ++ Consortium " , _email = " feelpp -devel@feelpp . org " )); auto mesh = unitSquare (); // define the mesh auto Vh = Pch <1 >( mesh ); // function space auto u = Vh -> element (); // element of function space auto v = Vh -> element (); // element of function space auto a = form2 ( _trial = Vh , _test = Vh ); // bilinear form //a = Ω ∇u • ∇v a = integrate ( _range = elements ( mesh ) , _expr = gradt ( u )* trans ( grad ( v )) ); auto l = form1 ( _test = Vh ); // linear form // l = Ω v l = integrate ( _range = elements ( mesh ) , _expr = id ( v )); // apply u = 0 on ∂Ω a += on ( _range = boundaryfaces ( mesh ) , _rhs =l , _element =u , _expr = constant (0.) ); // solve the equation a . solve ( _rhs =l , _solution = u );

  Figure 1. RB Framework

Figure 2 .

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

  pa ra m e t e r s p a c e _ p tr ty pe Dmu ; //D auto mu_min = Dmu -> element (); auto mu_max = Dmu -> element (); // fill mu_min and mu_max to define " corners " of D mu_min << ... ; mu_max << ... ; // associate " corners " to D Dmu -> setMin ( mu_min ); Dmu -> setMax ( mu_max ); parameter_type mu ; // µ ∈ D auto Pset = Dmu -> sampling (); int e im_sam pling_size = 1 000; Pset -> randomize ( eim_sampling_size ); // nonlinear expression σ(T ) auto sigma = ref ( mu (0) )/( 1+ ref ( mu (1) )*( idv ( T ) -T0 ) ); // eim object auto eim_sigma = eim ( _model = solve ( g(U , V ; µ) = 0 ) , _element =T , // unkown needed for the evaluation of σ(T ) _parameter = mu , // µ _expr = sigma , // σ(T ) _space =X h , _name = " eim_sigma " , _sampling = Pset );// then we can have access to β coefficients std :: vector < double > beta_sigma = eim_sigma -> beta ( mu );

Figure 4 .

 4 Figure 4. A Bitter insert consists in a stack of copper alloy plates. Insulators are introduced in between plates to create a helicoïdal current path.Small holes are drilled to enable the water flow throught the magnet. Some tie-rods are inserted to ensure a good electrical contact between each plate by applying a pre-compression.

FigureFigure 6 .

 6 Figure 5. CRB simulation on a sector of a Bitter magnet

Figure 7 .

 7 Figure 7. On the left : an example of a polyhelix insert; on the right a detailled view of the inner radially cooled helix with a zoom on the cooling channels and insulators.

Figure 8 .

 8 Figure 8. Temperature field computed with reduced basis method on a radial helix and on a sector, for inputs [σ 0 : 5.01e + 7; α : 3.48e -3; L : 2.89e -8; h : 8.15e + 4; T w : 295.2]. The value of parameter V D corresponds to a given current of 25 kA.

Acknowledgements.

Acknowledgements

The authors wish to thank Vincent Chabannes for fruitful discussions during the Cemracs 2012. Stphane Veys and Christophe Prud'homme acknowledge the financial support of the project ANR HAMM ANR-2010-COSI-009 and FNRAE-RB4FASTIM. Finally the authors wish to thank the Cemracs 2012 and its organizers.