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

Pierre Ladevèze

Elías Cueto
email: ecueto@unizar.es

Pierre Ladeveze
email: pierre.ladeveze@lmt.ens-cachan.fr

A Short Review in Model Order Reduction Based on Proper Generalized Decomposition

come

Introduction

Today many problems in science and engineering remain intractable, in spite of the impressive progresses attained in mechanical modeling, numerical analysis, discretization techniques and computer science during the last decade, because their numerical complexity, or the restrictions imposed by different requirements (real-time, for instance) make them unaffordable for nowadays technologies.

We can enumerate different challenging scenarios for efficient numerical simulations:

-The first one concerns models that are defined in high dimensional spaces, usually encountered in quantum chemistry and kinetic theory descriptions of complex fluids. Models defined in high dimensional spaces suffer the socalled curse of dimensionality. If one proceeds to the solution of a model defined in a space of dimension N by using a standard mesh based discretization technique, where M nodes are used for discretizing each space coordinate, the resulting number of nodes reaches the astronomical value of M N . With M = 1000 (a very coarse description in practice) and N = 30 (a very simple model) the numerical complexity results 10 90 . It is important to recall that 10 80 is the presumed number of elementary particles in the universe! Quantum mechanics and molecular modeling of complex fluids are not the only branches of science that suffer from the curse of dimensionality. Consider for example a chemical process involving so few molecules of the reacting species that the use of the continuous concept of concentration is not valid. This situation is often found in genetic processes such as expression of genes. The state of such a discrete system is given by a probability distribution for the number of individual molecules of each one of the coexisting species. The balance equation governing the evolution of the system, the so-called chemical master equation, is again defined in a high-dimensional space that precludes a direct solution by means of standard mesh-based techniques.

There are of course alternative methods to address these high-dimensional problems indirectly, stochastic simulations being one of the foremost approaches. Stochastic techniques have their own challenges, however. While variance reduction is always an issue, it is quite difficult within the stochastic framework to implement parametric or sensitivity studies that go beyond the brute force approach of computing a large number of expensive, individual simulations.

-A second category of problems involves multiscale problems not necessarily defined in high-dimensional spaces, but whose spectrum of characteristic times or lengths is so wide that standard incremental discretization techniques cannot be applied. In such time-multiscale problems, for instance, the time step is extremely small as a consequence of numerical stability requirements. Thus, simulations over the much larger time interval of interest, which typically requires the solution of a large linear algebraic system at each time step, simply become impossible. Multiscale models involving a wide range of characteristic times abound in many fields. Reaction-diffusion models of the degradation of plastic materials, for example, describe chemical reactions occurring within microseconds coupled to diffusion of chemical substances taking place over years. In processes involving microwaves, ultrasounds, etc. or materials exhibiting different relaxation times, the difficulty related to time integration is crucial. The same scenario is found in solid mechanics where the constitutive equations are strongly non linear and coupled, involving many scales and different characteristic times. -Other challenging problems are defined in degenerated geometrical domains. By this we mean that at least one of the characteristic dimensions of the domain is smaller by several orders of magnitude than the others. This is the case of bar, plate or shell-like domains typical of structures or materials processing applications. In simple situations, such problems are readily transformed into reduced, one or two-dimensional approximate theories (e.g. the classical elastic plate theory). When geometrical or material non-linearities are present, however, it is usually impossible to derive lower-dimensional models of sufficient validity. Standard mesh-based discretization methods then quickly become impractical, in view of the compulsory discretization of the small length scales that yield extremely fine meshes. -Many problems in process control, parametric modeling, inverse identification, and process or shape optimization, usually require, when approached with standard techniques, the direct computation of a very large number of solutions of the concerned model for particular values of the problem parameters. Consider for example the optimization of a process where optimal parameter values must be determined for process operating conditions (e.g. speed, position and temperature of heaters) and material properties (e.g. thermal and rheological properties of the materials). Clearly, it would be useful to be able to simulate this process at once for all values of these parameters within a prescribed interval, and then perform data mining within this rather general solution to identify optimal values. -Traditionally, Simulation-based Engineering Sciences-SBES-relied on the use of static data inputs to perform the simulations. These data could be parameters of the model(s) or boundary conditions. The word static is intended to mean here that these data could not be modified during the simulation. A new paradigm in the field of Applied Sciences and Engineering has emerged in the last decade. Dynamic Data-Driven Application Systems (DDDAS) constitute nowadays one of the most challenging applications of simulation-based Engineering Sciences. By DDDAS we mean a set of techniques that allow the linkage of simulation tools with measurement devices for real-time control of simulations. DDDAS entails the ability to dynamically incorporate additional data into an executing application, and in reverse, the ability of an application to dynamically steer the measurement process. In this context, real time simulators are needed in many applications. One of the most challenging situations is that of haptic surgery simulators, where forces acting on the surgical tool must be translated to the peripheral device at a rate of 500 Hz. Control, malfunctioning identification and reconfiguration of malfunctioning systems also need to run in real time. All these problems can be seen as typical examples of DDDAS.

-Augmented reality is another area in which efficient (fast and accurate) simulation is urgently needed. The idea is supplying in real time appropriate information to the reality perceived by the user. Augmented reality could be an excellent tool in many branches of science and engineering. -Light computing platforms are appealing alternatives to heavy computing platforms that in general are expensive and whose use requires technical knowledge. One can imagine that the off-line parametric solution of many models will make possible the on-line manipulation of those general solutions by using very light computing platforms, as for example smartphones or tablets.

While the previous list is by no means exhaustive, it includes a set of problems with apparent no relationship among them that can be, however, treated in an unified manner as will be shown in this paper. Their common ingredient is our lack of capabilities (or knowledge) to solve them numerically in a direct, traditional way. In order to obtain a solution, some kind of model order reduction is thus compulsory.

Many years ago, in the 80's, Pierre Ladeveze proposed a separated representation of the space and time coordinates

u(x, t) ≈ i=Q i=1 X i (x) • T i (t) (1)
for performing efficient solutions of complex non-linear thermo-mechanical models. The radial approximation (1) was one of the main blocks of the powerful non-incremental and non-linear LArge Time INcrement (LATIN) solver. The corpus of literature devoted to this technique is vast, as proved in the book by Ladeveze on the topic [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-new approaches and non-incremental methods of calculation[END_REF], but remained in the form of space-time separations for many years.

A more general separated representation (in this case in the conformation space of complex fluids descriptions) was more recently employed by A. Ammar and F. Chinesta [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] for approximating the solution of multi-dimensional partial differential equations. A. Nouy considered also such separated representations for solving stochastic equations were the deterministic coordinates and the stochastic ones were separated, very much like in the radial, space-time approximation. Proper Generalized Decomposition-PGD-is the common name recently coined for techniques using such separated representations. The nature of the problem, the involved coordinates and the constructor of such approximations can be very different, as will be seen.

The general form of the separated representation involved in the PGD reads:

u(x 1 , . . . , x N ) ≈ i=Q i=1 F 1 i (x 1 ) • • • F N i (x N ) (2)
where x i denote a scalar or vector coordinate defined in a domain

Ω i of moderate dimension Ω i ⊂ R d , with d ≤ 3, in general.

Review on PGD Foundations and Applications

Proper Generalized Decomposition techniques allow to circumvent most of the challenges mentioned in the introduction in a unified way. They also make it possible to solve models never until now solved, or speeding up the solution of others up to a limit not yet reached by other techniques. However, because the technique is still in its infancy, many applications (possibly the most challenging and impressive ones) are still in progress, while many and mathematical foundations remain to be fully understood. See [START_REF] Chinesta | Recent advances and new challenges in the use of the Proper Generalized Decomposition for solving multidimensional models[END_REF] for a recent review. In what follows we are describing the PGD state of the art.

Advanced Non-linear Solvers Making Use of Space-Time Separated Representations

Even if non-linear models can be solved by using standard linearization strategies in the context of separated representations, similarly to [START_REF] Ammar | Non-incremental strategies based on separated representations: applications in computational rheology[END_REF] or [113], more robust strategies exist for addressing efficiently history-dependent non-linearities, making use of the natural separation of the local and nonlinear thermo-mechanical behavior from the global and linear equilibrium in which a space-time separated representation constitutes an essential ingredient allowing spectacular computing time savings. Thus, space-time separated representations can be viewed as a part of the LArge Time INcrement method (LATIN) which has been introduced by P. Ladeveze in 1985 [START_REF] Ladevèze | New algorithms: mechanical framework and development[END_REF][START_REF] Ladevèze | On a family of algorithms for structural mechanics[END_REF] and was the core of two books [START_REF] Ladevèze | Mécanique non linéaire des structures[END_REF][START_REF] Ladevèze | Nonlinear computationnal structural mechanics-new approaches and non-incremental methods of calculation[END_REF]. These books relate the story of the PGD until their publication. Originally designed to deal with nonlinearities, such as plasticity or viscoplasticity, in time dependent problems [19-21, 27, 40-42, 59-61, 68, 72], during the past three decades, many theoretical works have been done to enhance the strategy and to deal with fundamental problems in computational mechanics. The strategy has been applied in the context of large displacements in [START_REF] Boucard | A nonincremental approach for large displacement problems[END_REF][START_REF] Ladevèze | Mécanique non linéaire des structures[END_REF]. A domain decomposition technique has been proposed in [START_REF] Blanzé | Une méthode pour l'étude d'assemblage de structures massives[END_REF][START_REF] Blanzé | A modular approach to structure assembly computationsapplication to contact problems[END_REF][START_REF] Champaney | Large scale applications on parallel computers of a mixed domain decomposition method[END_REF][START_REF] Ladevèze | A large time increment approach with domain decomposition technique for mechanical non linear problems[END_REF][START_REF] Ladevèze | Formulation et stratégies "parallèles" pour l'analyse non linéaire des structures[END_REF] and multiscale features in space have been introduced in [START_REF] Ladevèze | A 2-level and mixed domain decomposition approach for structural analysis[END_REF][START_REF] Ladevèze | A new micro-macro computational strategy for structural analysis[END_REF][START_REF] Ladevèze | A micro-macro and parallel computational strategy for high heterogeneous structures[END_REF]. These features have been extended to time in [74-77, 79, 80, 100] and the scalability of the strategy illustrated in [START_REF] Ladevèze | The latin multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Néron | Proper generalized decomposition for multiscale and multiphysics problems[END_REF]110]. A multiparametric strategy for taking into account the variability parameters or to be used in identification techniques has been introduced in [2, [START_REF] Boucard | A multiple solution method for non-linear structural mechanics[END_REF][START_REF] Boucard | Une application de la méthode latin au calcul multirésolution de structures non linéaires[END_REF][START_REF] Boucard | Application of the LATIN method to the calculation of response surfaces[END_REF]. The strategy has been applied to multiphysics problems in [START_REF] Dureisseix | A multi-time-scale strategy for multiphysics problems: application to poroelasticity[END_REF][START_REF] Néron | A computational strategy for poroelastic problems with a time interface between coupled physics[END_REF][START_REF] Néron | A computational strategy for thermo-poroelastic structures with a time-space interface coupling[END_REF]. Some techniques specific to the manipulation of fields approximated by the PGD have been proposed in [65]. The problem of the verification of model reduction methods within the PGD was addressed in [START_REF] Ladevèze | On the verification of model reduction methods based on the proper generalized decomposition[END_REF].

The strategy was also successfully applied to many studies, strongly related to industrial problems: the prediction of damage in composites [START_REF] Allix | A damage prediction method for composite structures[END_REF][START_REF] Aubard | Damage modeling at two scales for 4D carbon/carbon composites[END_REF][START_REF] Ladevèze | Debond modeling for multidirectional composites[END_REF][START_REF] Ladevèze | Multiscale modeling and computational strategies for composites[END_REF]118], the computation of assemblies [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF], the simulation of dynamic shocks [START_REF] Boucard | Macro-meso models for joints submitted to pyrotechnic shock[END_REF][START_REF] Ladevèze | A modular approach to 3-d impact computation with frictional contact[END_REF][START_REF] Lemoussu | A 3d shock computational strategy for real assembly and shock attenuator[END_REF], the simulation of porous media [START_REF] Dureisseix | A latin computational strategy for multiphysics problems: application to poroelasticity[END_REF][START_REF] Néron | Accounting for nonlinear aspects in multiphysics problems: application to poroelasticity[END_REF] and the virtual testing of joints for the prediction of damping in space launchers [START_REF] Caignot | Virtual testing for the prediction of damping in joints[END_REF].

Multidimensional Models

As discussed in the introduction, some models are inherently defined in a configurational or conformation space that is high-dimensional in nature. Its discretization leads to the well-known curse of dimensionality if traditional, meshbased techniques are applied. Separated representations, on the contrary, provide with an efficient means of dealing with these requirements. It is thus natural that the origin of PGD techniques, has one of its two legs in this framework, the other being the space-time separation of the LATIN method.

Separated representations were applied for solving the multidimensional Fokker-Planck equation describing complex fluids in the kinetic theory framework. In [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], the authors addressed the solution of the linear Fokker-Planck equation describing multi-bead-spring molecular models in the steady state and in homogeneous flows. In that case the distribution function describing the molecular conformation only depends on the molecule configurational coordinates. The solution procedure was extended to non-linear kinetic theory descriptions of more complex molecular descriptions in [START_REF] Mokdad | On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach[END_REF]. The transient solution was addressed in [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] in which the time was added as an extra-coordinate. Transient solutions making use of reduced bases in the context of an adaptive proper orthogonal decomposition [116] were considered in the case of low dimensional configuration spaces: the FENE model was addressed in [START_REF] Ammar | On the reduction of kinetic theory models related to finitely extensible dumbbells[END_REF] and liquid crystalline polymers in [START_REF] Ammar | Reduced numerical modeling of flows involving liquid-crystalline polymeres[END_REF]. In [START_REF] Leonenko | On the solution of the Fokker-Planck equation using a high-order reduced basis approximation[END_REF], the authors considered multi-bead-spring models but considered spectral approximation for representing all the functions involved in the finite sums decomposition. A deeper analysis of non-linear and transient models was considered in [START_REF] Ammar | Non-incremental strategies based on separated representations: applications in computational rheology[END_REF]. Complex fluid models were coupled with complex flows in [112] and [START_REF] Mokdad | A fully deterministic micro-macro simulation of complex flows involving reversible network fluid models[END_REF] opening very encouraging perspectives and claiming the necessity of defining efficient stabilizations. A first tentative of convective stabilization was proposed in [43]. Finally, in [START_REF] Chinesta | On the reduction of stochastic kinetic theory models of complex fluids[END_REF] the PGD was applied for solving the stochastic equation within the Brownian Configuration Field framework. The interested reader can refer to [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF] and the references therein for an exhaustive overview of PGD in computational rheology.

Multidimensional models encountered in the finer descriptions of matter (ranging from quantum chemistry to statistical mechanics descriptions) were revisited in [START_REF] Ammar | The nanometric and micrometric scales of the structure and mechanics of materials revisited: an introduction to the challenges of fully deterministic numerical descriptions[END_REF]. The multidimensional chemical master equation was efficiently solved in [START_REF] Chinesta | On the use of proper generalized decompositions for solving the multidimensional chemical master equation[END_REF]. The Langer's equation governing phase transitions was considered in [START_REF] Lamari | On the solution of the multidimensional Langerõs equation by using the proper generalized decomposition method for modeling phase transitions[END_REF]. Finally, models coming from financial mathematics were addressed in [START_REF] Falco | Algorithms and numerical methods for high dimensional financial market models[END_REF].

Separating the Physical Space

Other models are not inherently defined in a high-dimensional space, but can nevertheless be treated efficiently in a separated manner. Models defined in cubic (or hypercubic, but moderate-dimensional) domains suggest the following separated representation

u(x, y, z) ≈ i=Q i=1 X i (x) • Y i (y) • Z i (z) (3)
Thus, the 3D solution results in a sequence of 1D solutions for computing functions X i (x), Y i (y) and Z i (z). This kind of domains are very frequent when addressing homogenization problems in which the elementary volume element is a cube. The interested reader can refer to [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF]. A full decomposition was also efficiently applied for solving the Navier-Stokes equations in a cube in [START_REF] Dumon | Proper general decomposition (PGD) for the resolution of Navier-Stokes equations[END_REF]. Fully separated representations in complex, non-hypercubic domains can be equally performed by using an appropriate X-FEM like strategy [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF]. In that work the enforcement of nonhomogeneous Dirichlet boundary conditions was also analysed, giving rise to a generalized form of the PGD approach.

In the case of plates, shells or extruded geometries, one could consider the separated approximation [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF] also advantageously:

u(x, y, z) ≈ i=Q i=1 X i (x, y) • Z i (z) (4)
In the case of models defined in plate geometries (x, y) represent the in-plane coordinates and z the thickness direction. In the case of extruded profiles (x, y) represents the surface extruded in the z direction. This kind of representation makes possible fully 3D solutions with numerical complexities (and therefore computation times) characteristic of 2D solutions, without any simplifying kinematic assumption.

Parametric Models: A Route to Efficient Optimization, Inverse Identification and Real Time Simulation

Previous multidimensional models may seem too far from the usual computational mechanics practice. However, usual computational mechanics models could be enriched by a PGD treatment. Thus, adding some new coordinates to models initially non high-dimensional, could lead to new, never explored insights in the physics. Imagine for instance solving the heat (Fourier) equation and assume unknown the material thermal conductivity. This could happen because it has a stochastic nature, for instance, or simply because prior to solve the thermal model it should be measured. You have three possibilities: (i) wait for the conductivity to be measured (a conservative solution, but impractical in many engineering situations); (ii) solve the equation for many values of the conductivity (a sort of brute force approach) to get an overall idea of the behavior of the solution; or (iii) solve the heat equation only once for any value of the conductivity.

Obviously the third alternative is the most exciting one, but up to our knowledge, no result has been presented in the literature so far that could eventually allow to do so. To compute this general solution, however, it suffices to introduce the conductivity k as an extra-coordinate, playing the same role as the standard space and time coordinates, even if there are not derivatives concerning this extra-coordinate in the heat equation. A possible separated approximation could be:

u(x, t, k) ≈ i=Q i=1 X i (x) • T i (i) • K i (k) (5)
This procedure works very well in practice, and can be extended to introduce many other extra-coordinates: source term, boundary conditions, initial condition, etc. The price to pay is the increase of the dimensionality of the resulting model that now contains the standard physical coordinates (space and time) plus all the other extra-coordinates that we decided to introduce. However, the separated representation of the PGD allows to treat these models advantageously.

This kind of parametric modeling was addressed in [START_REF] Ammar | Solving parametric complex fluids models in rheometric flows[END_REF][START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF]114] where material parameters were introduced as extra-coordinates. In [START_REF] Lamari | Nonconventional numerical strategies in the advanced simulation of materials and processes[END_REF] and [START_REF] Lamari | Routes for efficient computational homogenization of non-linear materials using the proper generalized decomposition[END_REF] thermal conductivities, macroscopic temperature and its time evolution were introduced as extra-coordinates for computing linear and nonlinear homogenization. In [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF] the anisotropy direction of plies involved in a composite laminate were considered as extra-coordinates. By assuming a certain uncertainty in the real orientation of such plies, authors evaluate the envelope of the resulting distorted structures due to the thermomechanical coupling.

The futurist concept of "virtual chart for dimensioning" consists of a PGD-reduced model describing the physics for a family of structures. It allows engineers to proceed to complex designs in some seconds. This concept, fully compatible with the most recent design procedures, is being actively developed by the Ladeveze's group.

Real-Time Simulation, DDDAS and More

It is easy to understand that after performing this type of calculations, in which parameters are considered advantageously as new coordinates of the model, a posteriori inverse identification or optimization can be easily handled. This new PGD framework allows to perform this type of calculations very efficiently, because in fact all the result have been previously computed in the form of a separated, high-dimensional solution so that they constitute a simple post-processing of this general solution. Process optimization was considered in [START_REF] Ch | Methodological approach to efficient modeling and optimization of thermal processes taking place in a die: application to pultrusion[END_REF], for instance. Shape optimization was performed by considering all the geometrical parameters as extra-coordinates, leading to the model solution in any of the geometries generated by the parameters considered as extra-coordinates [START_REF] Leygue | A first step towards the use of proper general decomposition method for structural optimization[END_REF]. This strategy could be an alternative to the POD-based shape optimization considered in [117]. Inverse methods in the context of real-time simulations were addressed in [START_REF] Gonzalez | Proper generalized decomposition based dynamic datadriven inverse identification[END_REF] and were coupled with control strategies in [START_REF] Ch | Proper generalized decomposition based dynamic datadriven control of thermal processes[END_REF] as a first step towards DDDAS (dynamic data-driven application systems). Moreover, because the general parametric solution was precomputed off-line, it can be used on-line under real time constraints and using light computing platforms like smartphones [START_REF] Bognet | Advanced simulation of models defined in plate geometries: 3D solutions with 2D computational complexity[END_REF][START_REF] Ch | Proper generalized decomposition based dynamic datadriven control of thermal processes[END_REF], that constitutes a first step towards the use of this kind of information in augmented reality platforms.

As mentioned before, surgical simulators must operate at frequencies higher than 500 Hz. The use of model reduction seems to be an appealing alternative for reaching such performances. However, techniques based on the use of the POD, PODI, even combined with an asymptotic numerical methods to avoid the computation of the tangent matrix, exhibit serious difficulties to fulfil such requirements as discussed in [START_REF] Niroomandi | Real-time deformable models of non-linear tissues by model reduction techniques[END_REF][START_REF] Niroomandi | Model order reduction for hyperelastic materials[END_REF][START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced order models[END_REF][START_REF] Niroomandi | Real time simulation of surgery by reduced order modeling and X-FEM techniques[END_REF]. The use of parametric solutions in which the applied load and its point of application are considered as extra-coordinates open an unimaginable field of applications.

Uncertainty Quantification and Stochastic Parametric Analyses

The PGD method was introduced in [101] in the context of uncertainty quantification. In this context, the method was initially called Generalized Spectral Decomposition (GSD) and was considered to be a generalization of Karhunen-Loève decomposition for the a priori construction of separated representations of the solution of stochastic partial differential equations (SPDE). The problem that defines the reduced bases of functions in the decomposition of the solution has been interpreted as a pseudo eigenproblem. This interpretation has led to the development of dedicated algorithms inspired from solution techniques for classical eigenproblems [103]. In this context of stochastic problems, a major interest of PGD algorithms is that they decouple the solution of deterministic problems and algebraic stochastic equations, thus making the PGD a promising alternative to traditional methods for uncertainty propagation [105]. PGD has also been extended to stochastic nonlinear problems with two alternative strategies: the use of nonlinear iterative solvers with the PGD method as a solver for the succession of linearized stochastic problems [102], or the direct construction of a separated representation of the solution of nonlinear problems [104]. More recently, the PGD has been successfully applied to the solution of high dimensional stochastic parametric problems, with the introduction of suitable hierarchical tensor representations and associated algorithms [106].

PGD has also been exploited for the construction of robust reduced order models in structural solid mechanics [START_REF] Chevreuil | Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics[END_REF], thus providing an alternative to traditional modal superposition approaches. In [109], PGD method was coupled with fictitious domain methods in order to handle the solution of PDEs with geometrical uncertainties.

Numerical Analysis

From the point of view of the numerical analysis of separated representation and its associated constructors we can cite some recent works. There exist several approaches to the numerical analysis of PGD. They combine the existence of a best approximation and a greedy algorithm. The idea of using greedy algorithms to construct successive approximations was considered in [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] in the context of the numerical solution of multidimensional PDEs. One can note, however, that the strategy considered in [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] can be interpreted as a greedy algorithm for the error minimization only for elliptic symmetric problems since they can be rewritten in terms of minimization of an energy functional. In [START_REF] Ammar | On the convergence of a greedy rank-one update algorithm for a class of linear systems[END_REF] the analysis of the convergence of the greedy rank-one update algorithm for solving full-rank linear systems was addressed.

A more complete theoretical study was addressed in [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF] for symmetric elliptic problems by exploiting its connection to greedy algorithms from nonlinear approximation theory explored, for example, in [START_REF] Devore | Some remarks on greedy algorithms[END_REF]. Hence, the variational version of the algorithm, based on the minimization of a sequence of Dirichlet energies, was shown to converge. In [START_REF] Figueroa | Greedy approximation of highdimensional Ornstein-Uhlenbeck operators with unbounded drift[END_REF] the authors extend the convergence analysis of the pure greedy and orthogonal greedy algorithms considered in [START_REF] Bris | Results and questions on a nonlinear approximation approach for solving high-dimensional partial differential equations[END_REF] to the technically more complicated case where the Laplace operator is replaced by a high-dimensional Ornstein-Uhlenbeck operator with unbounded drift. This case bears many similarities with that appearing in Fokker-Planck equations arising in bead-spring chain type kinetic polymer models with finitely extensible non-linear elastic potentials. These are frequently posed on a high-dimensional Cartesian product configuration space

D = D 1 × • • • × D N contained in R Nd , where each set D i , i = 1, . . . , N, is a bounded open ball in R d , d = 2, 3.
General elliptic problems were considered in [START_REF] Falcó | A proper generalized decomposition for the solution of elliptic problems in abstract form by using a functional Eckart-Young approach[END_REF]. The extension of these results to non-linear symmetric coercive problems is the subject of the paper [START_REF] Cancès | Convergence of a greedy algorithm for high-dimensional convex nonlinear problems[END_REF] that also contains specific results when the algorithm is applied to finitedimensional problems. Current research tracks include for example the extension of the convergence analysis to nonsymmetric problems and the analysis of the rate of convergence for non-linear equations. In [START_REF] Falco | Proper generalized decomposition for nonlinear convex problems in tensor Banach spaces[END_REF], it is provided a mathematical analysis of a family of progressive and updated PGDs for a particular class of convex optimization problems in reflexive tensor Banach spaces. The above results have led to a generalization of the concept of Singular Value Decomposition (SVD). In particular, a constrained version of the SVD has been proposed in [108].

Miscellaneous

There is a variety of open questions when applying the PGD based discretization techniques. We discussed previously the question related to the convective stabilization [43]. The issue related to the stabilization of mixed formulation has not been addressed until now.

Other important question lies in the error estimation. A first attempt was considered in [START_REF] Ammar | An error estimator for separated representations of highly multidimensional models[END_REF]. Coupling standard mesh based discretization techniques with reduced bases (POD and PGD) is an issue of major interest in order to represent localized behaviors (e.g. discontinuities, singularities, boundary layers, . . .). Some coupling strategies were proposed in [START_REF] Ammar | Coupling finite elements and reduced approximation bases[END_REF] and [START_REF] Ammar | Coupling finite elements and proper generalized decompositions[END_REF].

Multi-scale and multi-physics non-linear coupled models involving different characteristic times were efficiently coupled in [START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF] where a globalization of local problems to speed-up the simulation in a separated representation framework has been proposed. Later, in [START_REF] Ammar | Proper generalized decomposition of time-multiscale models[END_REF] the time axis was transformed in a two-dimensional domain to account for two very different time scales. The coupling of other multiphysic models in the context of composites manufacturing processes was performed in [113]. A monolithic approach of coupled thermo-mechanical models was employed in [START_REF] Beringhier | Solution of strongly coupled multiphysics problems using space-time separated representations: application to thermoviscoelasticity[END_REF]. In [111] the solution of the electromagnetic equations was considered.

In [START_REF] Bonithon | Nonincremental boundary element discretization of parabolic models based on the use of proper generalized decompositions[END_REF] the PGD was introduced in the boundary element method framework for solving transient models where few works concerning the use of reduced bases exist [115]. Because the use of a space-time separated representation the time-dependent kernel is no more needed, and only the steady state kernel functions are needed. The spacetime separation implies significant computational time savings.

In [107], alternative definitions of PGD and dedicated algorithms are provided in the context of time-dependent PDEs. In particular, a new definition of PGD called Minimax PGD is introduced for handling non symmetric variational problems.

Discussion

As can be noticed from the pages before, the Proper Generalized Decomposition technique has been advantageously applied to a plethora of different problems and in different ways. Ranging from the original space-time separated representations to more recent applications in which models are efficiently cast into a high-dimensional framework to overcome real-time requirements by just making real-time postprocessing, the potential field of applications seems not to have been completely explored.

In fact, some crucial aspects at the very root of the method continue to be not fully understood. For instance, the relationship among POD/SVD and PGD techniques, has not been analysed to the whole extent. Optimality conditions of the separated representation, for instance, are intimately related to this question and have not been completely solved so far.

We believe, therefore, that there is a lot of room in the PGD field to perform research, both at the theoretical level and from the point of view of applications, whose range seems to be unlimited.

This work has been partially supported by the IUF-Institut Universitaire de France-, the french ANR COSINUS SIMDREAM project and the Spanish Ministry of Science and Innovation, through grants number CICYT-DPI2008-918 and DPI2011-27778-C02-01.

F. Chinesta ( )


David González

Elías Cueto

Francisco Chinesta

Pedro Díez

Antonio Huerta

Streamline upwind/Petrov-Galerkin-based stabilization of proper generalized decompositions for high-dimensional advection-diffusion equations

Keywords: convection-diffusion equation, SUPG, proper generalized decomposition, separated representations, finite sum decomposition, Fokker-Planck equation

This work is a first attempt to address efficient stabilizations of high dimensional advection-diffusion models encountered in computational physics. When addressing multidimensional models, the use of mesh-based discretization fails because the exponential increase of the number of degrees of freedom related to a multidimensional mesh or grid, and alternative discretization strategies are needed. Separated representations involved in the so-called proper generalized decomposition method are an efficient alternative as proven in our former works; however, the issue related to efficient stabilizations of multidimensional advectiondiffusion equations has never been addressed to our knowledge. Thus, this work is aimed at extending some well-experienced stabilization strategies widely used in the solution of 1D, 2D, or 3D advection-diffusion models to models defined in high-dimensional spaces, sometimes involving tens of coordinates.

INTRODUCTION

Despite the continuous growth of computer performance, some models remain nowadays intractable for current, state-of-the-art simulation techniques. Models defined in highly dimensional spaces, for instance, constitute a clear example of such problems. The interested reader can consult, for instance, [START_REF] Keunings | Micro-macro methods for the multiscale simulation of viscoelastic flow using molecular models of kinetic theory[END_REF] for an excellent survey of the state of the art on numerical simulation of micro-macro models for complex fluids. These micro-macro models are very often based on the Fokker-Planck equation, an advection-diffusion equation that describes the evolution of the probability distribution function of the phase-space variables, namely D Dt D @ @X ¹A º C 1 2 @ @X @ @X W D

where .x, X , t/ represents the probability distribution function of the phase-space variables X . Whereas A is an N d -dimensional deterministic vector responsible for the drift of these variables by the macroscopic flow, and D represents an N d N d -dimensional stochastic matrix responsible for Brownian effects on the flow. For Rouse models, for instance, N d can easily take values on the order of 50 or more [START_REF] Owens | Computational Rheology[END_REF].

Recently, the authors have presented and applied to different classes of problems an a priori model reduction technique based upon the approximation of the essential field by a finite sum of separable functions that are generated 'on the fly' by the method itself, see, for instance [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 representations[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[END_REF][START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | The nanometric and micrometric scales of the structure and mechanics of materials revisited: an introduction to the challenges of fully deterministic numerical descriptions[END_REF][START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] González | Recent advances on the use of separated representations[END_REF]. This method, coined as proper generalized decomposition (PGD, as opposed to a posteriori methods as the proper orthogonal decomposition [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Probability Theory[END_REF][START_REF] Lorenz | Empirical Orthogonal Functions and Statistical Weather Prediction[END_REF]), also allows for an easy treatment of problems defined in spaces of a high number of dimensions. Mesh-based techniques lead to the well-known curse of dimensionality, because the number of degrees of freedom grows exponentially with the number of dimensions. The use of separated representations allows to overcome this difficulty and has been used in many fields of Science and Engineering, such as in the framework of quantum chemistry. In particular, the Hartree-Fock (that involves a single product of functions) and post-Hartree-Fock approaches (as the multi-configurational self-consistent field that involves a finite number of sums) make use of a separated representation of the wave function [START_REF] Cancès | Computational quantum chemistry: a primer[END_REF]. In the context of Computational Mechanics, a similar decomposition was proposed, which was called radial approximation, and that was applied for separating the space and time coordinates in thermomechanical models [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF].

A detailed analysis of the convergence properties of different approaches of the PGD technique can be found at [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]. But this analysis has been restricted, for the moment, to symmetric operators. It remains unknown, to the best of the authors' knowledge, how the PGD behaves, for instance, when dealing with convection-diffusion equations.

It is well-known that standard finite element (Galerkin) methods do not work well for convectiondiffusion or convection-diffusion-reaction equations, because they lead to unstable, oscillating solutions [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF]. The first stabilization methods in finite elements including upwinding of convective terms that eventually lead to stable solutions were published in [START_REF] Christie | Finite element methods for second order differential equations with significant first derivatives[END_REF][START_REF] Heinrich | An 'upwind' finite element scheme for two-dimensional convective transport equation[END_REF].

Among the very numerous methods that have been proposed for the stabilization of convectiondiffusion equations, the streamline upwind/Petrov-Galerkin (SUPG) method [START_REF] Hughes | A multidimensional upwind scheme with no crosswind diffusion[END_REF] is one of the most extended. When a reaction term is important, subgrid scale (SGS) techniques have been advocated to reduce oscillations, see [START_REF] Hughes | The variational multiscale method-a paradigm for computational mechanics[END_REF]. An inherent difficulty of these methods is the choice of the stabilization parameter. In fact, algebraic or asymptotic analyses have been developed in one-dimensional problems. Optimal stabilization parameters in higher dimensions are not easily obtained.

This paper has neither the aim of addressing a state of the art on stabilization techniques for advection-diffusion equations and the consequent evaluation and choice of the optimal stabilization technique nor the proposal of new stabilization strategies. This work is a first attempt to address efficient stabilizations of high dimensional advection-diffusion models encountered in computational physics. When addressing multidimensional models, the use of mesh-based discretization fails because the exponential increase of the number of degrees of freedom related to a multidimensional mesh or grid and alternative discretization strategies are needed. Separated representations involved in the so-called PGD method is an efficient alternative as proven in our former works, some of them already cited; however, the issue related to efficient stabilizations of multidimensional advection-diffusion equations has never been addressed. Thus, this work has, as main aim, the extension of some well-experienced stabilization strategies widely used in the solution of 1D, 2D, or 3D advection-diffusion models to models defined in high dimensional spaces, sometimes involving tens of coordinates (the space and the time but also a number of configurational coordinates describing some extra model features). Thus, in what follows, we will focus on the extension of standard SUPG or subgrid scale stabilizations to solutions expressed in a separated approximation format.

Here, the use of finite sums of separable functions will be used to reduce a high-dimensional problem to a sequence of one-dimensional ones where optimal stabilization could be employed. The outline of the paper is as follows. In Section 2, we review the technique of separation of variables developed by the authors in the framework of convection-diffusion problems. In Section 3, we describe the framework of the class of convection-diffusion-reaction problems we deal with.We then provide two alternative formulations and discuss the advantages of each one in Sections 4 and 5. In Section 6, we show some examples illustrating the performance and capabilities of the technique so as to verify their stabilizing properties. These techniques are then applied in Section 7 to a 5-dimensional convection-diffusion problem.

THE PROPER GENERALIZED DECOMPOSITION

The use of a separated representation for solving complex models is not new. It has been employed in a variety of settings, including, as mentioned before, quantum mechanics in the Hartree-Fock approach or continuum mechanics, see the pioneering work by Ladeveze and coworkers [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF]. It is based on the establishment of an approximation of the unknown field, u, in the form:

u.x/ u n .x 1 , x 2 , : : : , x n sd / D n X iD1 F i 1 .x 1 / F i 2 .x 2 / : : : F i n sd .x n sd / , ( 1 )
with n sd as the number of spatial dimensions of the problem (n sd 1 within the scope of this work, although the results are general for any number of dimensions). The term proper generalized decomposition has been coined for methods that use this kind of separated representations, making use of an appropriate algorithm to determine the functions F i j . This a priori character of the technique is opposed to the a posteriori character of other model reduction techniques as proper orthogonal decomposition or Karhunen-Loeve method [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF][START_REF] Loève | Probability Theory[END_REF][START_REF] Ryckelynck | On the a priori model reduction: overview and recent developments[END_REF].

For the sake of simplicity in the description of the technique and without loss of generality, twodimensional problems are considered (they will be generalized later on, see Section 6); thus, the unknown field can be written as

u.x, y/ u n .x, y/ D n X iD1 F i .x/ G i .y/. ( 2 )
To construct this separated representation, an iterative, greedy algorithm is proposed. Assuming that the first n functions in the sum have been already computed, see [START_REF] Owens | Computational Rheology[END_REF], the method proceeds by finding the best functions R.x/ and S.y/ such that the updated representation given by

u nC1 .x, y/ D n X iD1 F i .x/ G i .y/ C R.x/ S.y/ D u n .x, y/ C R.x/ S.y/,
satisfies the weak formulation for a particular choice of the test functions. These test functions are w.x, y/ D R .x/ S.y/ C R.x/ S .y/.

(

) 3
Once these test functions are introduced in a weak form, such as [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF] or [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF], the resulting nonlinear problem of finding functions R and S is solved by using an iterative, fixed-point method. First, it is assumed that R is known and a new approximation for S is computed. Then assuming S is known, a new approximation for R is sought. This process is repeated until convergence.

For more details on the separated representation constructor, the interested reader can refer to [START_REF] Chinesta | Recent advances in the use of the proper generalized decomposition for solving multidimensional models[END_REF][START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] and the references therein.

PROBLEM SETTING

For simplicity and without loss of generality, we consider the steady-state convection-diffusionreaction equation. This equation is given by

a r u r . r u/ C u D s in R n sd , ( 4 )
with n sd the number of spatial dimensions-in the applications we are interested in much larger than one-and with boundary conditions

u D u D on D , ( 5 a
)

n r u D @u @n D t on N , ( 5b
)
where u is the scalar unknown field, a is the advective velocity, > 0 the diffusivity, assumed (without loss of generality) constant, the reaction term and s.x/ a volumetric source term. The function u D denotes the prescribed value of u on the Dirichlet portion of the boundary given by D , and t denotes the value of the normal diffusive flux on the Neumann boundary N .

The weak form of the problem defined in Equations ( 4) and ( 5)

is find u.x/ 2 S D ¹u 2 H 1 . /ju D u D on D º such that for all w 2 V D ¹w 2 H 1 . /jw D 0 on D º Z w.a r u/d C Z r w . r u/d C Z w ud D Z wsd C Z N wtd. (6)
This is usually expressed compactly as

a.w, u/ C c.aI w, u/ C .w, u/ D .w, s/ C .w, t/ N , ( 7 )
where the following definitions are employed:

a.w, u/ D Z r w . r u/d , c.aI w, u/ D Z w.a r u/d , .w, u/ D Z wud , .w, t/ N D Z N wtd.
The general form of a consistent stabilization technique is [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF] a

where P.w/ is some operator applied to the test functions, R.u/ D L.u/ s is the residual of the equation, and is the stabilization parameter. Note that the differential operator L is the one associated to the strong form, Equation (4), namely

L.u/ D a r u r . r u/ C u.
In the SUPG method, P.w/ D a r w, whereas in SGS, P.w/ D L .w/, where L is the adjoint operator, that is,

L .u/ D a r u r . r u/ C u.
Nodally exact results are obtained for linear elements in 1D for the convection-diffusion if

D h 2a  coth P e 1 P e à , ( 9 )
where P e is the mesh Péclet number, defined as P e D ah=2 ; h represents the mesh size parameter and a the modulus of the convective velocity. For convection-diffusion-reaction a fourth-order accurate formula reads

D h 2a 1 C 9 P e 2 C  h 2a à 2 ! 1=2 . ( 10
)
As noted previously, generalizations to higher dimensions are not readily available [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF].

The purpose of the remainder of this paper is to introduce a stabilized technique based on the use of the PGD introduced in [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 representations[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[END_REF] for high-dimensional problems. This kind of approximation makes use of a representation of the essential field in terms of a finite sum of separable functions.

In Sections 4 and 5, two alternative approaches based on PGD schemes are developed. A discussion with numerical examples follows. The first approach does the separation in infinite dimensional spaces, then when the finite element discretization is needed, the optimal stabilization parameter is chosen. The second one proceeds conversely, first the high-dimensional convection-diffusion equations is discretized, taking into account that stabilization is needed, and then separation following the PGD rationale is imposed.

STABILIZATION OF THE CONTINUOUS PGD SEPARATION

The PGD method is applied to the infinite dimensional problem defined by the weak form, Equation [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF]. In what follows, for the sake of simplicity, the reaction term of the equation is neglected, D 0, because it does not imply any special difficulty for the proposed methodology. Thus, the weak problem ( 7) becomes find R.x/ and S.x/ for all w or more precisely from (3) for all R .x/ and S .x/ in the proper infinite dimensional spaces, such that a.w, RS / C c.aI w, RS / D .w, s/ C .w, t/ N a.w, u n / c.aI w, u n /.

As noted earlier, the PGD method proceeds by a sort of alternating direction strategy, assuming iteratively that either R or S is known. For instance, in the previous equation, assuming that R is given and that S must be found, the weak form for all S reads a.RS , RS / C c.aI RS , RS / D .RS , s/ C .RS , t/ N a.RS , u n / c.aI RS , u n /, which is, in fact, a 1D convection-diffusion-reaction problem (recall that is standard to assume in PGD that D x ˝ y ), namely

S 0 , y S 0 y C S , c y S 0 y C S , Ä y S y D S , b y y C S p y ˇ ‡y Ny , (11)
where the following definitions are used:

y D Z x R 2 dx, c y D Z x a y R 2 dx, Ä y D Z x R 0 2 C a x RR 0 Á dx, b y D Z x  sR R 0 @u @x n R a r u n à dx C d dy ÄZ x R @u @y n dx , and
p y D Z ‡ x Ny t R dx Z x R @u @y n dx. (12)
It is important to remark that given the usual PGD assumption D x ˝ y , Neumann boundary conditions are also separable in the x and y problems. Here, because the unknown is S.y/, if N \@ y D N y ¤ ;, then N y must be separable as N y D ‡ y N y ˝ ‡ x N y . Thus, p y defined in (12) requires to integrate for a fixed y (corresponding to the Neumann boundary, denoted symbolically as ‡ y N y ) along the x direction, which is described as ‡ x N y . Note that proceeding identically with 'known' S and test function R S, an equivalent problem to determine R can be stated. This would induce a problem similar to Equations ( 11) and ( 12) but along the x direction.

In any case, Equation ( 11) can be identified as a convection-diffusion-reaction weak form, which is stabilized as follows:

where the operator applied to the test functions is P y .w/ D c y w 0 in SUPG and P y .w/ D c y w 0 C y w 0 0 Ä y w in SGS. Moreover, optimal values of the stabilization parameter can be employed substituting diffusion, y , convection, c y , and reaction, Ä y , in Equations ( 9) or [START_REF] Karhunen | Uber lineare methoden in der wahrscheinlichkeitsrechnung[END_REF].

For linear elements, the second derivatives in the stabilization terms are zero but recall that for high-order elements, consistent stabilization is important [START_REF] Huerta | Time accurate consistently stabilized mesh-free methods for convection-dominated problems[END_REF].

SEPARATION OF THE DISCRETE STABILIZED EQUATION

The other natural possibility arising from the PGD method is to perform a separated approximation of the already stabilized Equation ( 8). The method proceeds in a very similar way, but now, only one stabilization parameter must be determined (as in standard finite element techniques). The same stabilization parameter will be used in each one-dimensional problem. Moreover, although in each 1D direction is a different convection-diffusion-reaction problem, will be unique and defined a priori. Thus, in principle, for a convection-diffusion problem ( D 0) such as the one studied here, reaction will not be accounted for in the stabilization.

In the stabilized weak problem [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF], once the PGD method proceeds assuming iteratively that either R or S is known, becomes, when for instance, R is given and that S must be found, To compare with the previous formulation, the 1D problem in y for SUPG is presented. However, such comparison is only possible if further assumptions are imposed, namely, that diffusion and convection are independent of x and y, which is not always the case. Under these assumptions, the previous equation can be rewritten as

S 0 , y S 0 y C S , c y S y C S , Ä y S y C X e Z e y a y S 0 c y S 0 y S 00 C Q Ä y S dy C X e Z e y a x S O c y S 0 O y S 00 C O Ä y S dy D S , b y y C S p y ˇ ‡y Ny X e Z e a x R 0 S C a y RS 0 .a r u n r . r u n / s/ d , ( 14
)
where the definitions of ( 12) are used and the following ones are needed

Q Ä y D Z x RR 00 C a x RR 0 dx, O y D Z x RR 0 dx, O c y D Z x a y RR 0 dx, and O Ä y D Z R 0 R 00 C a x R 0 2 Á dx.
x This clearly evidences that both formulations are not, even in the simplified case, nearly equivalent. Moreover, it is obvious that the stabilization of this 1D problem is not optimal. Nevertheless, as shown in the examples, stabilized results are obtained.

VERIFYING THE STABILIZATION PROCEDURES

Before addressing the solution of a high-dimensional advection-diffusion equation, we are verifying the previous extensions of standard stabilization towards the separated representation format work. For this, we are focusing in some well-illustrated 2D benchmark problems, to which a reference solution can be easily computed to conclude on the performances of the considered stabilizations within the PGD framework. After proving the ability of these procedures for addressing the stabilization within the PGD framework, we will address a problem defined in a high-dimensional space to which no reference solution is available.

A first example

In this section, we consider an example of convection-diffusion skew to the mesh, proposed in [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF], among other references, see Figure 1.

The flow is unidirectional, with jjajj D 1, with the convective velocity not aligned to any of the axes, forming 30°with the x-direction. Thus, the results will be entirely comparable with those in reference [START_REF] Donea | Finite Element Methods for Flow Problems[END_REF]. Following this same reference, the diffusivity coefficient is taken to be 5 10 6 , corresponding to a mesh Péclet number of 10 4 . The inlet boundary conditions are discontinuous and of two different types at the outlet: We compare the results obtained with formulations 1 (introduced in section 4) and 2 (section 5) explained before in the framework of proper generalized decompositions. In Figure 2, results for the Dirichlet problem solved with PGD and standard Galerkin approach are shown. No stable results were obtained for a mesh of 2000 finite elements along each spatial direction.

Downwind homogeneous natural boundary conditions Downwind homogeneous essential boundary conditions

To compare with existing techniques, the problem has been solved by SUPG-stabilized standard finite elements. A refined mesh of 100 100 finite elements has been employed to obtain a reference solution to compare with. In addition, and for the purpose of comparison, a 10 10 mesh was employed. These results are depicted in Figures 3 and4 for Dirichlet and Neumann boundary conditions, respectively.

Results obtained with the PGD formulation 1 and 2 are depicted in Figure 5 for the Dirichlet problem and Figure 6 for the Neumann problem. As it can be noticed, results are practically indistinguishable for this problem. These results were obtained with the product of 85 functions for version 1 and 20 functions for version 2 along each direction. These functions are depicted in Figure 7 for the Dirichlet problem and formulations 1 and 2. The great number of very similar functions obtained by the PGD method and version 1 can be noticed. Note that in the aforementioned Figure 7, the approximating functions have been normalized. Although 85 versus 20 functions seems to be a great difference for both approaches to the problem, by looking at the convergence plots in Figure 8, it can be noticed that both formulations, for a fixed number of summands, give very similar levels of accuracy. It seems, however, that the functions incorporated by the second approach are closer to be mutually orthogonal in comparison with those of version 1. Namely, their cosines, computed as their scalar product divided by the product of their norms, is lower for version 2 than for version 1.

To clearly compare the results, a cut along x D 0.5 has been carried out. These are plotted in Figures 9 and10 Taking the 100 100 SUPG-stabilized finite element mesh results as a reference, we can compute the norm of the 'error' with respect to this reference. These values are reported in Tables I and II for the Dirichlet and Neumann boundary conditions, respectively.

As can be noticed, the proposed formulation (both, in fact) provides a somewhat more steep slope in the jump while there appears to be slightly high oscillations. To better understand the approximation provided by the proposed method, we include the next example.

A second example

We consider the problem cos 3

Á @u @x C sin 3 Á @u @y D 10 4 Â @ 2 u @x 2 C @ 2 u @y 2 Ã C 1,
solved on the unit square, with boundary conditions (see Figure 11)

u D 1 on x D 0, u D 1 on x D 1, u D 0 on y D 0, u D 0 on y D 1.
In the absence of analytical solution, a FEM-SUPG solution on a mesh composed by 100 100 linear elements is employed as a reference, see Figure 12.

Results are compared with the reference solution and that of a 10 10 standard linear finite elements with SUPG stabilization. As a reference value, for the formulation 2, 20 functions along each spatial direction were employed, see Figure 13. Similar values are obtained for formulation 1.

Results are better analyzed by comparing a cut of the obtained results at abscissa x D 0.8. These are shown in Figure 14 below. If we compare the obtained errors with respect to the reference solution, we obtain the results summarized in Table III.

A third example

This time we solved a problem in the unit square with kak D 5, oriented at 30°with the horizontal, and we took D 1, see Figure 15. Homogeneous Dirichlet boundary conditions were assumed on the whole boundary. Again, a 100 100 finite element mesh, with standard SUPG stabilization, was considered as a reference solution. In this example, the differences between the proposed method and the equivalent in finite elements are more evident. Table IV resumes these differences. The reference solution and the obtained approximations for formulations 1 and 2 are depicted in Figure 16. Note that the PGD approaches to the problem given in this example has a remarkable difference in accuracy with the standard, SUPG-stabilized finite element approach. Note also that the first version of section 4 gives more accurate results in the L 2 -norm than that in section 5.

Also noticeable is the number of functions needed for attaining such level of error. The algorithm coined as 'version 2' needed only one pair of functions to give the error shown in Table IV. For one pair of functions, version 1 of the method gave an L 2 -norm error of 0.01645, still one order of magnitude less than the stabilized FEM.

A HIGH DIMENSIONAL EXAMPLE

To fully justify not only the behavior of the proposed technique, which has been thoroughly addressed in the previous sections, but also its true advantages for high-dimensional problems, we consider here a problem defined in R 5 .

The problem is defined in the hypercube .x, y, ´, t , s/ 2 OE0, 1 5 , with a D .1, 1, 1, 1, 1/. All boundary hyperplanes were subjected to Dirichlet boundary conditions of the same type: u D 0 at x, ´, t , s D 0, u D 1 at x, ´, t , s D 1, whereas u D 1 at y D 0 and u D 0 at y D 1. A finite element mesh of ten elements per edge of the hypercube was employed. The restricted number of elements along each direction is kept here for the purpose of highlighting the stabilization characteristics of the proposed techniques rather than its ability to cope with the curse of dimensionality, which has been addressed in previous works, see [START_REF] González | Recent advances on the use of separated representations[END_REF], for instance. Obviously, the representation of results in a 5D domain is not easy. We have chosen to represent the results at the plane .x, y/ located at ´, t , s D 0.5. Obviously, the solution without stabilization is highly oscillating, as expected, see Figure 17.

On the contrary, results for the two different stabilization techniques presented in Sections 4 and 5, respectively, are shown in Figure 18. Functions employed to approximate the solution along each spatial direction (for the formulation in Section 4 in this case) are depicted in Figure 19. These same results can be now compared with those obtained by employing a mesh of 100 elements along each edge of the hypercube. This would make a total of 100 5 D 10 10 elements if traditional finite element methods would have been employed. Obviously, this is out of reach for a standard personal computer (this is a practical consequence of the curse of dimensionality mentioned in the introduction). However, the result in Figure 20 has been obtained by means of PGD approximations in less than 2 minutes in a laptop running MATLAB. Mathworks, Inc. 3 Apple Hill Drive Natick, MA 01760-2098 UNITED STATES.

CONCLUSIONS

The PGD techniques allow for an efficient means of a priori model reduction. While constructing alternatively a reduced approximation to the problem by means of products of separable functions, PGD techniques have been devised as an efficient method to deal with the so-called curse of dimensionality, that is, the exponential growth of the number of degrees of freedom for mesh-based discretization techniques applied to problems defined in highly-dimensional spaces.

Fokker-Planck equations, for instance, constitute a typical example of such difficulty, because they are usually defined in highly dimensional spaces. Fokker-Planck equations are a class of conservation (convection-diffusion) equations for the probability density function associated with the time evolution of the micro-state variables in multi-scale models of complex fluids.

In this work, we have addressed precisely an analysis of the behavior of PGD techniques in the context of convection-diffusion equations. Two main possibilities arose. The first one, more appealing at a first sight, consists in establishing a separated representation of the infinite-dimensional convection-diffusion equation. This gave rise to a sequence of one-dimensional convectiondiffusion-reaction problems that could eventually be stabilized by SUPG techniques or whatever state-of-the-art stabilization technique suitable for this class of problems. SUPG techniques present the advantage of exact stabilization (at the nodes) for one-dimensional convection-reaction problems solved with linear finite elements. Thus, PGD techniques seemed to be specially well-suited for this problem, not only providing a reduced modeling technique for the problem but also attaining exact stabilization in any number of dimensions. This approach gave good results, as discussed in the text.

The second approach consisted in applying the PGD technique to the already stabilized problem and showed a similar degree of accuracy as formulation 1. The price to pay, however, is the multidimensional (and therefore not exact) stabilizing term to be added to the equation before constructing the approximation. This technique, however, greatly exploited its dimensional reduction capabilities and showed great promise for its extension to highly-dimensional problems. For the examples included in this paper, only 20 terms were enough for approximating the problem with a resolution that would make necessary to employ 10 10 degrees of freedom with standard finite elements.

  .w, u/ C c.aI w, u/ C .w, u/ C X e Z e P.w/ R.u/d D .w, s/ C .w, t/ N ,

S 0 ,

 0 y S 0 y C S , c y S y C S , Ä y S y C X e Z e y P y .S / y c y S 0 y S 0 0 C Ä y S Á dy D S , b y y C S p y ˇ ‡y Ny C X e Z e y P y .S / y b y dy,

  a.RS , RS / C c.aI RS , RS / C X e Z e P.RS / .a r .RS / r . r .RS /// d D .RS , s/ C .RS , t/ N a.RS , u n / c.aI RS , u n / X e Z e P.RS / .a r u n r . r u n / s/ d .

Figure 1 .

 1 Figure 1. Statement of the problem used for the validation of the proposed technique.

Figure 2 .

 2 Figure 2. Solution obtained by applying the separated representation approach to standard, non stabilized Galerkin procedures.

Figure 3 .

 3 Figure 3. Dirichlet problem: FEM-SUPG result with a (a) 10 10 mesh and (b) 100 100 mesh.

Figure 4 .Figure 5 .

 45 Figure 4. Neumann problem: FEM-SUPG result with a (a) 10 10 mesh and (b) 100 100 mesh.

Figure 6 .

 6 Figure 6. Neumann problem: results of formulations 1 and 2.

Figure 7 .

 7 Figure 7. Functions employed in the approximation of problem 1 with Dirichlet boundary conditions, (a) along x-coordinate and (b) y-coordinate. Formulations 1 (top row) and 2 (bottom row) are depicted.

  for the Dirichlet and Neumann problems, respectively.

Figure 8 .

 8 Figure 8. Convergence for the Dirichlet problem and formulations 1 and 2. Error in L 2 -norm versus number of functions incorporated to the approximation.

Figure 9 .Figure 10 .

 910 Figure 9. Comparison between standard, finite element stabilized (SUPG) approximations and the proposed PGD formulations for the Dirichlet problem.

Figure 12 .

 12 Figure 12. Reference solution for the problem in section 6.2, obtained by FEM-SUPG and a 100 100 linear finite element mesh.

Figure 13 .Figure 14 .

 1314 Figure 13. Comparison between the obtained results by (a) FEM-SUPG (10 10 elements), (b) formulation 1, and (c) formulation 2.

Figure 16 .Figure 17 .

 1617 Figure 16. Comparison between the obtained results by (a) FEM-SUPG (100 100 elements, reference solution), (b) FEM-SUPG (10 10 elements) (c) formulation 1, and (d) formulation 2.

Figure 18 .Figure 19 .Figure 20 .

 181920 Figure 18. Comparison between the obtained results by (a) formulation 1 and (b) formulation 2.

Table I .

 I Comparison of 'errors' with respect to the reference solution of a 100 100 SUPG-stabilized finite element mesh. Dirichlet problem.

	Formulation	1	2	FEM-SUPG
	jjejj L 2 jjejj H 1	0.249436 1.249441	0.240834 1.242036	0.273473 1.273479
	Table II. Comparison of 'errors' with respect to the reference solu-
	tion of a 100 100 SUPG-stabilized finite element mesh. Neumann
			problem.
	Formulation	1	2	FEM-SUPG
	jjejj L 2 jjejj H 1	0.251757 1.251761	0.241130 1.242258	0.259762 1.259767

Figure 11

. Geometry of the problem analyzed in section 6.2 .

Table III .

 III Obtained 'errors' with respect to the 100 100 SUPG-stabilized finite element solution.

	Formulation	1	2	FEM-SUPG
	jjejj L 2 jjejj H 1	0.016660 1.016674	0.0101530 0.999251	0.085855 1.09428
		Figure 15. Geometry of problem 3.

Table IV .

 IV Obtained 'errors' with respect to the 100 100 SUPG-stabilized finite element solution for problem 3.

	Formulation	1	2	FEM-SUPG
	jjejj L 2 jjejj H 1	0.009163 1.159664	0.009871 1.009154	0.1596 1.159664

ACKNOWLEDGEMENTS

This work is supported by the Spanish Ministry of Economy and Innovation CICYT-DPI2011-27778-C02-01/02.