Amine Ammar
email: amine.ammar@ujf-grenoble.fr

Etienne Prulière

Julien Férec

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Elías Cueto
email: ecueto@unizar.es

Etienne Pruliere

Elias Cueto

Coupling Finite Elements and Reduced Approximation Bases

Keywords: model reduction, proper orthogonal decomposition, separated representations, finite elements. MOTS-CLÉS : réduction de modèles, POD, séparation de variables, éléments finis

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Most engineering systems can be modeled by a continuous model usually expressed as a system of linear or non-linear coupled partial differential equations describing the different conservation balances (momentum, energy, mass and chemically reacting substances). From a practical point of view, the determination of the problem solution, that is, the knowledge of the different fields characterizing the physical system at any point and time (velocity, pressure, temperature, concentrations…) is not possible in real systems due to the complexity of models, geometries and boundary conditions. For this reason the solution is searched only at some points and at some times, from which it could be interpolated to any other point and time. Techniques leading to this kind of representation are known as discretization techniques. There are today numerous discretization techniques (finite elements, finite volumes, finite differences, meshless techniques…). The optimal technique to be applied depends on the model and on the domain geometry. The development of numerical analysis and the computation availabilities make possible today the resolution of complex systems involving millions of unknowns related to the discrete model. However, the complexity of the models is also increasing exponentially, and today engineers are not only interested in solving models, but also in solving these modes very fast and, if possible, with a great accuracy.

In the context of control, optimization, inverse analysis, and in general computation in real time, it is clear that numerous problems must be solved, and for this reason the question related to the computation time becomes crucial. The question is very simple: is it possible to perform very fast and accurate simulations? Different answers have been given to this question depending on the scientific community to which this question is addressed. For specialists in computational science an answer to this question requires the improvement of computational resources, high performance computing and the use of parallel computing platforms. For some specialists in numerical analysis the challenge is in the fast resolution of linear systems via the use of preconditioners or multigrid techniques among many others. For others the idea is to adapt the cloud of nodes (points where the solution is computed) in order to avoid excessive number of unknowns. However, today all these approaches allow to alleviate the computation efforts but the fast and accurate computation remains a real challenge. This paper presents another different approach, allowing impressive reduction in the number of degrees of freedom without compromising the computed solution accuracy, based on the use of reduced approximation bases. The idea is very simple. Consider a domain where a certain model is defined and the associated cloud of nodes able to represent by interpolation the solution everywhere. In general the number of unknowns scales with the number of nodes, and for this reason even if the solution is evolving in time smoothly all the nodes are used for describing it at each time step. In the reduced modeling that we are describing in this paper the numerical algorithm is able to extract the optimal information describing the evolution of the solution in the entire time simulation interval. Thus, the evolution of the solution can be expressed as a linear combination of a reduced number of functions (defining the reduced approximation basis), and then the size of the resulting linear problems are very low, and the CPU time savings of some orders of magnitude can be attained (Ammar et al., 2006).

The extraction of this relevant information is a well known topic that can be performed by applying the proper orthogonal decomposition -POD -, also known as Karhunen-Loève decomposition -KLD - (Karhunen, 1946;[START_REF] Loève | Probability theory[END_REF] that will be summarized in the next section. This kind of approach has been widely used for weather forecast purposes [START_REF] Lorenz | Empirical orthogonal functions and statistical weather prediction[END_REF], turbulence [START_REF] Sirovich | Turbulence and the dynamics of coherent structures part I: Coherent structures[END_REF][START_REF] Holmes | Lowdimensional models of coherent structures in turbulence[END_REF], solid mechanics [START_REF] Krysl | Dimensional model reduction in non-linear finite element dynamics of solids and structures[END_REF] but also in the context of chemical engineering for control purposes (see the Park's works, e.g. [START_REF] Park | The Use of the Karhunen-Loève decomposition for the modelling of distributed parameter systems[END_REF].

Usual reduced models perform the simulation of some similar problems or the desired one in a short time interval. From these solutions the POD decomposition can be performed, allowing extracting the most relevant functions describing the solution evolution. Now, it is assumed that the solution of a "similar" problem can be expressed using this reduced approximation basis, allowing to a significant reduction on the discrete problem size and then to significant CPU time savings. However, in general the question related to the accuracy of the computed solutions is usually ignored. An original approach combining the model reduction and the control of the solution accuracy was proposed by [START_REF] Ryckelynck | A Priori Hyperreduction Method: an adaptive approach[END_REF], and applied later in different domains [START_REF] Ryckelynck | On the 'a priori' model reduction: Overview and recent developments[END_REF]Ammar, Ryckelynck et al., 2006;[START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF].

In the adaptive technique initially proposed by [START_REF] Ryckelynck | A Priori Hyperreduction Method: an adaptive approach[END_REF], the so called "a priori" model reduction, the solution of a transient model is expressed, under certain regularity requirements, from a finite sum decomposition, that is, as the sum of some products, each one involving a function of space and a function of time. For this purpose, one must define a procedure to define the functions of space, and the ones depending on time are computed by enforcing the verification of the dynamical system. The main issue of this kind of technique is precisely the procedure to buildup the optimal reduced approximation basis. Ryckelynck proposed the use of the residual or some Krylov's subspaces generated by the residual to enrich the approximation basis, and after convergence, the just enriched basis is reduced by applying the proper orthogonal decomposition (the algorithm will be summarized later). However, some times (e.g. the wave equation with a moving load) the enrichment procedure based on the use of the Krylov's subspaces doesn't run.

Another possibility to define a separated representation of a solution consists of building-up each term in the finite sum decomposition from the Galerkin form of the governing equation by using an alternating direction scheme, looking for the fixed point of the iteration solver. This strategy was originally introduced by Ladeveze (see Ladeveze 1999 and the references therein) that called it "radial approximation". This strategy was later and independently used for the solution of stochastic models [START_REF] Nouy | Recent developments in spectral stochastic methods for the solution of stochastic partial differential equations[END_REF] and by [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: Transient simulation using space-time separated representation[END_REF] for addressing highly multidimensional models. The main advantage of this technique with respect to the one based on the use of the adaptive POD is that both functions (the one defined in space and the one defined on time) are built-up by the algorithm without any major assumption in a transparent way for the user. Moreover, the first numerical experiments and the first numerical analysis results proved its superior robustness. Thus, both techniques just presented can be used for defining reduced approximation bases able to compute fast and accurately the solution of transient models. However, some times the solutions of general models could involve localization, singularities, discontinuities (fixed or evolving)… and then the coupling between reduced approximation bases (for describing the smooth part of the solution) and standard finite elements for describing the localized behavior seems a very appealing strategy in computational mechanics. This work is a first attempt towards this ambitious objective.

Even if significant computing time savings are expected when using this coupled strategy, this paper only focus in the analysis of the coupling strategy and the evaluation of the reduction of the number of degrees of freedom involved in the resulting discrete model. A deep analysis of the computing time savings will be addressed in an oncoming paper.

A priori model reduction based on the use of the POD

We assume that the evolution of a certain field depending on the physical space . This is equivalent to obtaining functions x maximizing

2 1 1 2 1 p P i N p i i p i i N i i u x x x [1]
which leads to:

1 1 1 1 p P j N i N i N p p i i j j i i p i j i u u x x x x x x [2]
, which can be rewritten in the form

1 1 1 1 j N p P i N i N p p i j j i i i i j p i u u x x x x x x [3]
Defining the vectors a such that its i-component is i a x , Equation [3] takes the following matrix form

; T T k k [4]
where the two points correlation matrix k is given by 1 1

p P p P T p p p p ij i j p p u u k x x k u u [5]
which is symmetric and positive definite. If we define the matrix Q containing the discrete field history:

1 2 1 1 1 1 2 2 2 2
1 2

P P P N N N u u u u u u u u u Q [6]
is easy to verify that the matrix k in Equation [5] results

T k Q Q [7]

Reduced modeling

If the evolution of a certain field is known , , 1 , , , 1 , ,

p p i i u t u i N p P x [8]
which is coming from some direct simulations or from experimental measures, then matrices Q and k can be computed and the eigenvalue problem given by Equation [4] solved. The solution of Equation [4] results in N couples of eigenvalueeigenvector. However, in a large number of models involving smooth time evolutions of the solution, the magnitude of the eigenvalues decreases very fast, evidencing that the solution evolution can be represented as a linear combination of a reduced number of functions (the ones related to the highest eigenvalues).

In our numerical applications we consider the eigenvalues ordered n are selected, because their associated eigenvectors are enough to represent accurately the entire solution evolution. In a large variety of models n N and moreover n only depends on the regularity of the solution evolution, and neither on the dimension of the physical space (1D, 2D or 3D) nor on the size of the model N. We verified in the numerical examples addressed in this work but also in many other addressed in former works that by considering a threshold value of 8 10 the error between the reduced solution and fully finite element solution (L 2 norm) was in the order of 4 10 . If the threshold is reduced to 2 10 for example, error of about 0.1 (L 2 norm) were noticed in the examples addressed in this paper. Even if we cannot conclude, some numerical experiments seem to indicate that the error scales with the square root of the threshold value.

The reduced approximation basis consists on the n eigenvectors 1 , , n , allowing to define the basis transformation matrix B :

1 2 , , , n B [9]
whose size is N n . Thus, the vector containing the field nodal values u can be expressed by: 1 ( ) ( )

n i i i a t t u B a [10]
Now, if we consider the linear system of equations resulting from the discretization of a partial differential equation (PDE) in the form

Enriching the reduced approximation basis

The just described strategy allows for very fast computation of large size models. For example one could solve the full model using some standard discretization technique (finite differences, finite elements…) for a small time interval and then define matrix Q and k allowing compute the reduced approximation basis transformation B that leads to the reduced solution procedure illustrated by Equation [13].

However, it is not guaranteed that this reduced basis that has been built from the solution evolving in a short time interval remains accurate for describing the solution in the entire simulation interval. Moreover, during the simulation material properties, boundary conditions… could change, and sometimes these changes are significant. Thus, in general the level of confidence expected for this reduced model solution is decreasing as the time (out the interval that served to build the reduced approximation basis) increases.

In this manner, if one would compute reduced model solutions and keep the confidence on the computed solution, a check of accuracy must be applied and, if needed, an enrichment strategy must be performed in order to adapt the reduced approximation basis for capturing the new events present in the solution evolution which can not be described accurately from the original reduced approximation basis. Using this strategy one can attain computing time savings of some orders of magnitude (Ammar et al., 2006).

For this purpose, [START_REF] Ryckelynck | A Priori Hyperreduction Method: an adaptive approach[END_REF] proposed to start with a low order approximation basis, using some simple functions (e.g. the initial condition in transient problems) or using the eigenvectors of a "similar" problem previously solved or the ones coming from a full simulation within a short time interval. Now, we compute S iterations of the evolution problem using the reduced model [13] without changing the approximation basis. After these S iterations, the complete discrete system [12] is constructed, and the residual R evaluated:

1 1 S S S S R = Au -f = ABa -f [14]
If the norm of the residual is small enough, R , with a threshold value small enough, we can continue for other S iterations using the same approximation basis. On the contrary, if the residual norm is too large, R , we need to enrich the approximation basis and compute again the last S iterations. This enrichment is built using some Krylov's subspaces, in our case the three first subspaces: 2 B B,R,AR,A R .

In the numerical example here addressed we considered 3 Krylov's subspaces (this choice is not essential) and start by putting S=1000. If an enrichment is needed after the S iterations, so we put S=S/2. If it is not the case, we put: S=2S.

One could expect that the enrichment process is increasing continuously the size of the reduced approximation basis, but in fact, after reaching the convergence, a POD decomposition is performed on the whole past time interval in order to extract the significant information. The reduced basis is then updated accordingly.

Obviously, because this control and enrichment stages we can ensure the quality of the computed solution.

We described in this section the main ideas of the adaptation strategy that we are using. However, in order to be efficient from a computational point of view, the numerical algorithm that we are using is a bit more complex. Interested readers can refer to [START_REF] Ryckelynck | On the 'a priori' model reduction: Overview and recent developments[END_REF] and the references therein for a more detailed and valuable description.

Towards a finite elements-reduced approximation bases coupling

Generic form of the model

In what follows we assume the physical model described by a generic partial differential equation where N is the vector containing the shape functions and matrix incorporates integration by parts and stabilizations of the advective terms eventually involved by the operator .

With this notation the discrete model results:

* * * 0, T T u u u u u u u 0 [17]

Defining the approximation coupling

This section deals with models involving localized singularities or discontinuities, being the solution far from these regions quite smooth. When these regions evolve in time, additional difficulties appear whose treatment is out of the scope of the present paper. Thus, one expects that efficient simulations could be achieved by coupling reduced descriptions in the smooth regions and a fully finite element description within the regions involving singularities or discontinuities.

Discontinuities can be efficiently captured in the extended finite elements framework, however the application of the techniques previously described in the context of the extended finite elements for describing evolving discontinuities introduces serious difficulties because the number of approximation functions needed becomes too large. For these reasons the simplest and most efficient solution consists of dividing the domain in two parts, and then applying different approximations in each of these resulting regions. The one that concerns the finite element approximation must cover the domain in which the singularity or discontinuity is evolving. Obviously, if the singularity (or discontinuity) explores the whole domain the technique that we are describing doesn't work, but in the most of engineering models these singularities remain confined in regions in general much smaller than the whole domain. Another difficulty related to the application of the coupling that we are describing is the necessity of an "a priori" knowledge about the region in which the evolving singularity (or discontinuity) will be confined.

Despite all the criticisms that we just addressed to the coupling strategy that we are describing in this section, numerous models accept this simple coupling, as we illustrate in the next section.

Let 2 be the small region of interest where localized singularities exist and 1 the remaining part of the domain (where the solution is assumed varying smoothly):

1 2

1 2

[18]

Each node of the discretization belongs to a single subdomain. The matrix form of the discrete model [17] is now rewritten, by introducing a decomposition of the unknown vector u in 1 u and 2 u , containing the degrees of freedom in the domains 1 and 2 respectively, as:

11 12 1 11 12 1 21 22 2 21 22 2 u u 0 u u 0 [19]
On the domain 1 the "a priori" model reduction described in the previous section could be applied. A set of reduced basis function are involved in the matrix B. For this purpose we introduce the following notation: Thus, at least in the region 1 the computing cost has been significantly reduced. The residual calculation is performed over the whole domain and for the enrichment of the basis, only the values of the residual on 1 are retained.

Because the use of a finite element interpolation of both the unknown function in region 2 (fully finite element description) and of the functions involved in the reduced approximation basis, the approximation continuity in the whole domain is ensured as when one uses a fully finite element description in the entire domain.

From an algebraic point of view, the band size increases as can be noticed in Equation [21]. However, as the reduced model contains very few degrees of freedom its impact on the linear solver is negligible.

Numerical examples

A simple heat transfer problem

To illustrate the just described technique we first consider a simple 2D heat transfer problem. This problem concerns the evolution of the temperature field ( , ) T t

x governed by the heat equation:

( ), , T T h t I t x x [22]
where the diffusivity is assumed isotropic and constant in the whole domain and time (it was fixed to a unit value in the simulations here reported),

( 2, 2) and 0,1 I . The heat source ( ) h x was taken in the form 9 9 9 9

( , ) max( , 0) max( , 0) min( , 0) min( , 0)

h x y x y x y [23]
in order to induce high gradient of the temperature field at the boundaries neighborhood.

Finally, the initial and boundary conditions are given by 2 2

( , , 0) cos ( ) cos ( ) ( , ) 0

T x y t x y T t x [24]
This problem was solved by employing a reduced approximation in the whole domain and also by applying the coupling technique described in the previous section. Figure 1 illustrates the solution computed by using the first strategy. The solution was approximated with only six couples of functions according to:

6 1 ( , ) ( ) ( ) i i i i T t F a t x x [25]
whose discrete counterpart writes:

6 1 ( ) ( ) i i i i t a t T F [26]
Figure 1 depicts the four most representative modes 1 4 ( ), , ( )

F F x
x , the evolution of all the reduced coordinates ( ), 1, , 6 i a t i as well as the resulting temperature field at 1 t . There are no appreciable differences between the finite element solution and the reduced one (the norm of the error is 4 2 10 ).

The solution employing the coupling method was performed by decomposing the physical domain en two regions, as depicted in Figure 2. In that figure the triangular elements located in the white region have nodes taking part of the standard finite element approximation and others that are concerned by the reduced approximation. The solution in the green region is performed by applying the model reduction previously described. To represent this solution only six modes were needed. The four most representative modes, the evolution of the reduced coordinates and the solution at time t = 1 are depicted in Figure 3. Finally Figure 4 depicts the finite element solution at six different times. By representing together the solutions at any time computed in the green and red regions, the solution at that time is defined in the whole domain.

Accounting for localized behaviour in kinetic theory models

The description of the microstructural conformation in complex fluids needs too many degrees of freedom. In this kind of models one must describe the molecular conformation that in the simplest descriptions (that correspond to the most widely used) is carried out by defining the molecule end-to-end vector x . Thus, a molecule is finally modelled as two beads connected by a spring. When this spring is Hookean it allows infinite extensions, however in order to limit the maximum molecule stretching different non-linear spring stiffness were introduced. One of the most popular is the so-called FENE model, in which the maximum extension is denoted by b and the spring stiffness writes:

2 1 ( ) 1 f b x x [27]
The microstructure state is then defined by introducing a probability distribution function (p.d.f.) ,t

x giving the fraction of molecules that at time t has the end-toend vector given by x . Obviously, because the finite extension of such molecules 0, B b x , that is, the end-to-end vector belongs to the ball of radius b .

The dimensionless conservation balance of the p.d.f.

1 1 ( ) 2 2 f v x x x x x [29] (a) (c) (e) (b) (d) (f) t [30]
and v is a given velocity gradient of the fluid assumed homogeneous in space. In what follows we are solving the just defined model by assuming a given flow whose kinematics is described from its gradient of velocity 0 6 0 0 v [33] being the parameter b=10. For the calculation of the stress related to the presence of molecules we are computing the virial stress (also known within the rheology community as the Kramer's rule) according to:

( ) ( ) ( , ) t f t d x x x x x [34]
Figure 5 depicts the domain decomposition. In the larger region a reduced modelling is performed, whereas in both small regions a fully finite element discretization is applied. We perform this decomposition because the solution behaviour is known in this case. However, when nothing on the solution is known an error indicator should be integrated in order to define the regions needing for a fully finite element description. The definition of appropriate error indicators constitutes a work in progress.

The computed solution making use of the reduced approximation basis can be described from a finite sums decomposition involving 14 terms:

14 1 1 ( , ) ( ) ( ) i i i i i t F G t x x [35]
where the space and time functions were normalized (by dividing by its norm). The four most representative space modes are depicted in Figure 6(a)-(d), where all the functions ( ) i G t are also represented. Finally the solution at time t=5 is reconstructed and shown in Figure 6(f). Figure 7 depicts the finite element solution in the remaining regions at different times. As in the example described in the previous section, by merging the solutions in both regions at different times, the solution in the whole domain is easily defined. Figure 8 depicts the steady solution in the whole domain and Figure 9 shows the evolution of the extra-stress due to the presence of molecules computed from [34].

In the model considered in this section the mesh contains 951 nodes. Obviously this number corresponds with the number of degrees of freedom if one is using a standard finite element model. However, 777 nodes are located in the region 1 that where reduced to 14 degrees of freedom involved in the reduced approximation. The remaining degrees of freedom (951-777=174) consist of the nodal values at nodes located in region 2 . Thus, the size of the final discrete systems scales with the number of nodes located in region 2 because the ones related to the reduced approximation (of the order of 10) are negligible with respect to the former ones.

Conclusions

In this paper we presented a first attempt for coupling reduced and finite element approximation bases. The reduced approximation basis allows significant reduction of the number of degrees of freedom involved in the resulting discrete model in the regions in which the solution evolves smoothly. On the other hand, in the regions where the solution exhibits high gradients or even discontinuities a standard finite element approximation seems better. Coupling of both descriptions is straightforward. Obviously no reduction is expected in the finite element description regions, however, in the remaining one, a reduced number of functions, and then of degrees of freedom. At present the user must specify the domain decomposition, but one can imagine that this task could be easily carried out automatically. The ongoing work consists of coupling separated and FE representations.

N

  . The n eigenvalues belonging to the interval 1 n with

  solution at the previous time step, now, introducing Equation [10] it results the final system of equations is of low order, i.e. the dimensions of T B A B are n n , with n N , and the dimensions of both a and Equation [13] can be also derived introducing the approximation [10] into the partial differential equation Galerkin form.

.

  and and are two differential operators. In what follows, and without loss in generality, we are assuming that t

Figure 1 .Figure 2 .

 12 Figure 1. Solution computed using a fully reduced approximation basis: (a)… (d) most representative modes; (e) reduced coordinates evolution; (f) solution at t = 1

Figure 3 .

 3 Figure 3. Solution in 1 : (a)… (d) most representative modes; (e) reduced coordinates evolution; (f) solution at t = 1

Figure 4 .

 4 Figure 4. Fully finite element solution in 2

Figure 5 .Figure 6 .

 56 Figure 5. Domain decomposition for the solution of the kinetic theory model

Figure 7 .Figure 8 .

 78 Figure 7. Solution in 2 : finite element representation