David Gonzalez

Amine Ammar

Francisco Chinesta

Elıas Cueto

Recent advances on the use of separated representations

Keywords: curse of dimensionality, separation of variables, singular value decomposition, essential boundary conditions

Separated representations based on finite sum decompositions constitute an appealing strategy for reducing the computer resources and the calculation costs by reducing drastically the number of degrees of freedom that the functional approximations involve (the number of degrees of freedom scale linearly with the dimension of the space in which the model is defined instead of the exponential growing characteristic of mesh-based discretization strategies). In our knowledge the use of separated representations is the only possibility for circumventing the terrific curse of dimensionality related to some highly multidimensional models involving hundreds of dimensions, as we proved in some of our former works. Its application is not restricted to multidimensional models, obviously separated representation can also be applied in standard 2D or 3D models, allowing for high resolution computations. Because its early life numerous issues persist, many of them attracting the curiosity of many research groups within the computational mechanics community. In this paper we are focusing in two issues never until now addressed: (i) the imposition of non-homogenous essential boundary conditions and (ii) the consideration of complex geometries.

INTRODUCTION

Some models encountered in science and engineering are sometimes defined in multidimensional spaces (as the ones involved in quantum mechanics or kinetic theory descriptions of materials, including complex fluids) that exhibit the terrific curse of dimensionality when usual mesh-based discretization techniques are applied.

In those models the difficulty is quite natural and their solution needs for new strategies. One possibility lies in the use of sparse grids [START_REF] Bungartz | Sparse grids[END_REF] but its use is restricted to models defined in spaces of moderate dimensions. Another technique able to circumvent, or at least alleviate, the curse of dimensionality consists of using a separated representation of the unknown field as we proposed in our former works [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 representations[END_REF] and then applied in numerous contexts: (i) quantum chemistry [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]; (ii) Brownian dynamics [START_REF] Chinesta | On the reduction of stochastic kinetic theory models of complex fluids[END_REF]; (iii) kinetic theory description of polymers solutions and melts [START_REF] Mokdad | On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach[END_REF]; (iv) kinetic theory descriptions of rods suspensions [START_REF] Pruliere | Multiscale modelling of flows involving short fibersuspensions[END_REF]; among many others.

Basically, the separated representation of a generic function u(x 1 , . . . , x D ) (also known as finite sums decomposition) writes:

u(x 1 , . . . , x D ) ≈ i=N i=1 F i 1 (x 1 )ו••× F i D (x D ) (1)
This kind of representation is not new, it was widely employed in the last decades 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 MCSCF that involves a finite number of sums) made use of a separated representation of the wavefunction [START_REF] Cancès | Computational Quantum Chemistry: A Primer, Handbook of Numerical Analysis[END_REF]. In the context of computational mechanics a similar decomposition was proposed, that 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].

This kind of approximation only needs a technique able to construct, in a way completely transparent for the user, the separated representation of the unknown field involved in a partial differential equation (PDE).

The technique that we proposed for computing the different functions involved in Equation (1) consists of an alternating directions linearization strategy that we summarize here. For the sake of clarity, and without any loss of generality, we restrict our discussion to the D-dimensional Poisson's equation

u = -f (x 1 , x 2 , . . . , x D ) ( 2
)
where u is a scalar function of (x 1 , x 2 , . . . , x D ). Problem (2) is defined in the domain (x 1 , x 2 , . . . , x D ) ∈ = (-L , +L) D with vanishing essential boundary conditions. The problem solution can be written in the form

u(x 1 , x 2 , . . . , x D ) = ∞ j=1 j D k=1 F k j (x k ) (3)
where F k j is the jth basis function, with unit norm, which only depends on the kth coordinate.

It is well known that the solution of numerous problems can be accurately approximated using a finite (sometimes very reduced) number (N ) of approximation functions, i.e.:

u(x 1 , x 2 , . . . , x D ) ≈ N j=1 j D k=1 F k j (x k ) (4)
The previous expression implies the same number of approximation functions in each dimension, but each one of these functions could be expressed in a discrete form using different number of parameters (nodes of the 1D grids). Now, an appropriate numerical procedure is needed for computing the coefficients j as well as the N approximations functions in each dimension.

The proposed numerical scheme consists of an iteration procedure that solves at each iteration n the following three steps:

Step 1: Projection of the solution in a discrete basis. If we assume the functions F k j (∀ j ∈ [1, . . . , n]; ∀k ∈[1, . . . , D]) known (verifying the boundary conditions), the coefficients j can be computed by introducing the approximation of u into the Galerkin variational formulation associated with Equation (2):

∇u * •∇u d = u * f d (5)
Introducing the approximation of u and u

* u(x 1 , x 2 , . . . , x D ) = n j=1 j D k=1 F k j (x k ) (6)
and

u * (x 1 , x 2 , . . . , x D ) = n j=1 * j D k=1 F k j (x k ) (7)
we have

∇ n j=1 * j D k=1 F k j (x k ) •∇ n j=1 j D k=1 F k j (x k ) d = n j=1 * j D k=1 F k j (x k ) f d (8)
Now, we assume that f (x 1 , . . . , x D ) can be written in the form

f (x 1 , . . . , x D ) ≈ m h=1 D k=1 f kh (x k ) (9)
Equation ( 8) involves integrals of a product of D functions each one defined in a different coordinate. Let D k=1 g k (x k ) be one of these functions to be integrated. The integral over can be performed by integrating each function in its definition interval and then multiplying the D computed integrals according to

D k=1 g k (x k ) d = D k=1 L -L g k (x k ) dx k ( 10
)
which makes possible the numerical integration in highly dimensional spaces. Now, due to the arbitrariness of the coefficients * j , Equation (8) allows to compute the n-approximation coefficients j , solving the resulting linear system of size n ×n. This problem is linear and moreover rarely exceeds the order of tens of degrees of freedom. Thus, even if the resulting coefficient matrix is densely populated, the time required for its solution is negligible with respect to the one required for performing the approximation basis enrichment (step 3).

Step 2: Checking convergence. From the solution of u at iteration n given by Equation ( 6) we compute the residual Re related to Equation (2):

Re = ( u + f (x 1 , . . . , x D )) 2 u ( 11
)
If Re< (epsilon is a small enough parameter) the iteration process stops, yielding the solution u(x 1 , . . . , x D ) given by Equation [START_REF] Mokdad | On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach[END_REF]. Otherwise, the iteration procedure continues.

The integral in Equation ( 11) can be written as the product of one-dimensional integrals by performing a separated representation of the square of the residual.

Step 3: Enrichment of the approximation basis. From the coefficients j just computed the approximation basis can be enriched by adding the new function D k=1 F k(n+1) (x k ). For this purpose we solve the non-linear Galerkin variational formulation related to Equation ( 2)

∇u * •∇u d = u * f d ( 12
)
using the approximation of u given by:

u(x 1 , x 2 , . . . , x D ) = n j=1 j D k=1 F k j (x k )+ D k=1 R k (x k ) (13)
The weighting function can be expressed as:

u * (x 1 , x 2 , . . . , x D )=R * 1 (x 1 )×R 2 (x 2 )ו••× R D (x D )+ • • • +R 1 (x 1 )×R 2 (x 2 )ו••× R * D (x D ) (14)
This leads to a non-linear variational problem, whose solution allows to compute the D functions R k (x k ). Functions F k(n+1) (x k ) are finally obtained by normalizing, after convergence of the nonlinear solver, the functions R 1 , R 2 , . . . , R D .

To solve this problem we introduce a discretization of those functions R k (x k ). Each one of these functions is approximated using a 1D finite element (FE) description. If we assume than p k nodes are used to construct the interpolation of function R k (x k ) in the interval [-L , L], then the size of the resulting discrete non-linear problem is k=D k=1 p k . The price to pay for avoiding a whole mesh in the multidimensional domain is the solution of a non-linear problem. However, even in high dimensions the size of the non-linear problems remains moderate and no particular difficulties have been found in its solution up to hundreds dimensions. Concerning the computation time, even when the non-linear solver converges quickly, this step consumes the main part of the global computing time.

Different non-linear solvers have been analyzed: fixed-point, Newton or one based on an alternating directions scheme. In this work the last strategy was retained. Thus, in the enrichment step, function R s+1 1 (x 1 ) is updated by assuming known all the others functions (given at the previous iteration of the non-linear solver R s

2 (x 2 ), . . . , R s D (x D )). Then, functions R s+1 1 (x 1 ), R s 3 (x 3 ), . . . , R s D (x D ) are assumed known for updating function R s+1 2 (x 2 )
, and so on until updating the last function

R s+1 D (x D ). Now the convergence is checked by calculating i=D i=1 R s+1 i (x i )-R s i (x i ) 2 .
If this norm is small enough we can define the functions F k(n+1) (x k ) by normalizing the functions R 1 , R 2 , . . . , R D and come back to step 1. On the contrary, if this norm is not small enough, a new iteration of the non-linear solver should be performed by updating functions R s+2 i (x i ), i = 1, . . . , D and then checking again the convergence. Despite its simplicity, our experience proves that this strategy is in fact very robust.

We must recall that the technique that we proposed in the papers just referred, is not a universal strategy able to solve any kind of multidimensional PDE. Thus, the efficient application of the technique that we just described requires the separability of all the fields involved in the model. Obviously, this separability is not always possible because some functions need a tremendous number of sums. On the other hand, even when the field is separable (one could perform this separation by invoking, for example, the singular value decomposition (SVD) or the multidimensional SVD) the finite sums decomposition of general multidimensional functions is not realistic because the amount of memory needed for storing the discrete form of such functions before applying the multidimensional SVD.

In many physical models (see for example [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 | 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]), a fully separation (consisting of a sum of products of one-dimensional functions) could not be envisaged from a practical point of view. Thus, a better approximation lies in writing

u(x 1 , . . . , x d ) ≈ i=N i=1 F i 1 (x 1 )ו••× F i D (x d ) ( 15
)
where the different functions taking part in the finite sums decomposition are defined in spaces of moderate dimensions, that is x i ∈ i ⊂ R q i , where usually q i = 1, 2 or 3.

When the model is separable, as it is the case in many physical models coming from computational science and engineering, the unknown field involved in the model can be expressed as a tensor product of approximations defined in lower-dimensional domains, eventually one-dimensional. However in that case, two major difficulties appear: (i) the first one is related to the enforcement of non-homogeneous essential boundary conditions; and (ii) the representation of geometrically complex domains other than cubes or hyper-cubes. This paper concerns these two difficulties that never until now, up to our knowledge, have been addressed. For the sake of clarity (and mainly for the ease of representation of the results), and without loss of generality, in what follows we restrict our analysis to 2D or 3D models.

The simplest way to enforce non-homogeneous essential boundary conditions (i.e. u = u g in * ) consists in defining a function û(x), x ∈ ⊂ R D , with the requested regularity, verifying the boundary conditions, i.e. û(x ∈ * ) = u g (x ∈ * ). Now, one could perform the separation of û, by invoking the SVD or the multidimensional SVD, leading to

û(x 1 , . . . , x D ) ≈ i=M i=1 G i 1 (x 1 )ו••×G i D (x D ) ( 16
)
that allows defining the approximation of the unknown field as

u(x 1 , . . . , x D ) ≈ i=M i=1 G i 1 (x 1 )ו••×G i D (x D )+ i=N i=1 F i 1 (x 1 )ו••× F i D (x D ) (17)
where the second part of the right member

i=N i=1 F i 1 (x 1 )ו••× F i D (x D ) is subjected to homogenous boundary conditions.
Finally, in some cases, even when the separability is assured, the geometrical complexity of the domain needs for special treatments. Different possibilities could be envisaged. This paper explores a very simple strategy based on the use of R-functions.

Sections 2 and 3 focus in the treatment of non-homogeneous essential boundary conditions and complex geometries, respectively, within the separated representation framework. Finally, Section 4 presents some numerical examples proving the efficiency of the strategies proposed.

IMPOSITION OF NON-HOMOGENEOUS ESSENTIAL BOUNDARY CONDITIONS

Consider, for simplicity and without loss of generality, the following problem:

-u = f (x 1 , x 2 , . . . , x D ) in = (0, L) D (18)
subjected to the boundary conditions:

u = u g = 0 on ≡ * (19)
Let us assume that we are able to find a function , continuous in , such that -∈ L 2 ( ) verifying Equation (19) [START_REF] Rektorys | Variational Methods in Mathematics[END_REF]. Then, the solution of the problem given by Equations ( 18) and ( 19) can be obtained straightforwardly by

u = + z (20)
where we thus face a problem in the z variable

-z = f + in (21) z = 0 on (22)
solvable by the method of separation of variables before presented.

Construction of the function

We thus take an extension of the essential boundary conditions, u g , , within the domain . The solution will be constructed in the form

u ≈ + Q j=1 j D k=1 F k j (x k ) (23)
thus achieving the exact imposition of essential boundary conditions. In the last equation, F k j (x k ) is a function of the approximation space of the solution with homogeneous boundary conditions. In practice, the function u g can be defined at different portions of the boundary, through different functions i . Assuming the domain where the model is defined given by the hypercube (0, L) D , essential boundary conditions must be enforced on all the (D -1)-dimensional faces of the hypercube: the two limit points of a 1D segment, the four faces of a 2D square, the six faces of a 3D cube, and so on. Thus, the number of faces of a D-dimensional hypercube results P = 2D on which different essential boundary conditions i are assumed prescribed.

u( 1 ) = 1 u( 2 ) = 2 . . . u( P ) = P ( 24
)
where i represent the different (D -1)-dimensional faces.

Each one of these functions i can be grouped to form a single, global, function utilizing, for instance, the transfinite interpolation method proposed in [START_REF] Rvachev | R-functions in boundary value problems in mechanics[END_REF]. Each portion of the essential boundary * i can be described by an implicit function w i (x 1 , x 2 , . . . , x i-1 , x i+1 , . . . , x D ) = 0, where w i is a function vanishing on * i and positive elsewhere (the distance function, for instance, can be considered as a suitable w i function). Thus,

= P i=1 i P j=1; j =i w j P i=1 P j=1; j =i w j (25)
constitutes a generalization of Lagrangian interpolation. The resulting function interpolates each i along each w i . Following this procedure, the weak form of the problem will be: Find z ∈ H 1 ( ) such that for every z * ∈ H 1 0 ( )

∇z * •∇z d = z * f d + ∇z * •∇ d (26)
holds.

Remark 1

The proposed form for the function , given by Equation ( 25) is singular at some points of the boundary. Provided that numerical integration has been performed in this work through standard, first order, Gauss quadrature rules, no integration point will be placed in these singular positions.

There exist, however, some alternative expressions for the function that posses singular points at different positions. The interested reader can consult their form in reference [START_REF] Rvachev | On completeness of RFM solution structures[END_REF] and references therein.

Finally, the obtained function , which is generally not product or sum separable, according to [START_REF] Messiter | Product and sum separable functions[END_REF], is separated approximately by employing the SVD technique. To this end, we employ singular values up to a given precision, which is normally fixed to 10 -5 .

PROBLEMS DEFINED OVER GENERAL DOMAINS

The proposed solution strategy is based on the method of the characteristic function [START_REF] Babuska | Survey of meshless and generalized finite element methods: a unified approach[END_REF]. Let us consider a function : → R, with a hypercube (usually the smallest one) containing the arbitrarily-shaped domain , and with continuous derivatives such that:

(x) > 0, x ∈ (x) = 0, x ∈ * (x) = 0, x / ∈ |∇ (x)| >0, x ∈ * (27)
The solution space will then be

S h = {u | u = v, v ∈ H 1 ( )}⊂ H 1 0 ( ).
The approximate solution to the problem given by Equations ( 18)-(19) will be u h ∈ S h , where u h is the solution of the problem:

Find u ∈ S h such that for every w ∈ H 1 0 ( )

∇w •∇u d = w f d ( 28
)
and S h = S h , with S h the finite-dimensional trial space, that depends on the particular method employed.

A study of the error and convergence associated with this method can be found in [START_REF] Babuska | Survey of meshless and generalized finite element methods: a unified approach[END_REF] and references therein, whereas the demonstration of the completeness of the method is accomplished in [START_REF] Rvachev | On completeness of RFM solution structures[END_REF].

The final form of the solution in the context of a separated representation will then be:

u ≈ • Q j=1 j D j=1 F k j (x k ) (29)
which will verify straightforwardly the homogeneous essential boundary conditions enforced along the boundary.

Remark 2

If the problem under consideration possesses non-homogeneous essential boundary conditions and is defined in a domain other than an hypercube at the same time, the obvious choice for the approximation combines the techniques developed in Sections 2 and 3 above:

u ≈ + • Q j=1 j D j=1 F k j (x k ) (30)

Construction of the function

It is clear that this approach is not very useful unless a method to easily construct such function in a general case is established. In this case we employ R-functions [START_REF] Rvachev | R-functions in boundary value problems in mechanics[END_REF].

An R-function is a real-valued function whose sign is completely defined by the sign of its arguments. Such functions encode boolean operations that help to construct the combinations of simple, basis, functions. For instance, consider the following functions [START_REF] Rvachev | R-functions in boundary value problems in mechanics[END_REF], that behave like the logical operators and and or:

x ∧ y ≡ (x + y -x 2 + y 2 ) (31) x ∨ y ≡ (x + y + x 2 + y 2 ) (32)
Consider, for instance, the domain shown in Figure 1. The domain is defined by means of a set of six inequalities of the general form These inequalities can be combined to render a definition of the whole domain as:

f 1 = y y 1 f 2 = x x 2 f 3 = y y 3 . . . ( 33
) f 1 f 2 f 3 f 4 f 5 f 6
= ( f 1 ∨ f 2 )∧ f 3 ∧ f 4 ∧ f 5 ∧ f 6 ( 34
)
An example of the performance of this technique is shown in Figure 2, where the essential boundary of the domain is defined (Figure 2 The question now lies in establishing a separated representation for the function . By employing standard SVD we can employ up to m singular values of the decomposition,

11 ו••× 1D , . . . , m1 ו••× m D
, necessary to obtain a prescribed precision. Thus, the final expression for the approximate solution of the problem will be

u ≈ [ 11 (x 1 )ו••× 1D (x D )+•••+ m1 (x 1 )ו••× m D (x D )]× N j=1 j D j=1 F k j (x k ) (35)
The difficulty in this case is that a large number of terms in the SVD decomposition of the function would lead to a complex expression for the approximate solution with a large number of terms in the sum. Each term would be integrated in its own 1D space through standard Gauss quadratures. Thus, finally, the above difficulty can be easily circumvented. which is represented in Figure 4. This function cannot be written as a closed sums decomposition according to [START_REF] Messiter | Product and sum separable functions[END_REF]. Applying SVD, we arrive to the one-dimensional functions depicted in Figure 5, that approximate the function to the chosen level of precision, 10 -5 in this case. The error generated in this separation, by employing the first six singular values of the decomposition, is depicted in Figure 6.

The method before described was applied in meshes composed by 15, 20, 30, 50, 80 and 130 nodes along each spatial direction. The solution was compared with one approximate solution obtained by standard FEs in a mesh composed by 40 000 nodes. The obtained solution, that will be considered as the reference one, is depicted in Figure 7.

The solution strategy can be resumed as follows. For each discretization, ranging from 15 nodes to 130 nodes at each spatial direction, we have employed the same SVD decomposition of function , constructed to approximate the function to an error of 10 -5 . The separated representation strategy presented in Sections 1 and 2 has been employed until reaching a small enough residual. The obtained error for each discretization along each spatial direction is depicted in Figure 8. The final solution for the separation of variables technique is depicted in Figure 9. The homogeneous part of the solution, i.e. the solution for the problem in z variable, is shown in Figure 10.

Extension of the proposed technique to higher-dimensional problems

The main issue associated with the straightforward extension of the proposed technique to higherdimensional problems lies in obtaining the appropriate counterpart of the SVD decomposition in high dimensions. To this end, we have employed the technique implemented by Kolda and Bader [START_REF] Kolda | Tensor decompositions and applications[END_REF], in conjunction with the MATLAB tensor toolbox [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF]. In general, we seek for decompositions of our function , whose nodal values are stored in an order n (n = 3 in the example here considered, for simplicity of graphical representation) tensor A, in the form

A = R× 1 U× 2 V× 3 W× 4 . . . ( 41
)
where R is a diagonal tensor and U, V, W, . . . are the components or factors of the decomposition.

Columns of the factors are orthonormal in standard, two-dimensional SVD, but this is not possible, in general, for tensors of order three or higher. The products involved in the previous decomposition are defined as follows. Mode-n multiplication, denoted by × n , means to multiply the matrix times the tensor in the nth mode. In other words, the n-mode product of a tensor are given by

C ∈ R I 1 ×I 2 ו••×I N by a matrix U ∈ R J n ×I n is an I 1 × I 2 ×. . . I n-1 × J n × I n+1 ו••× I N -
(C× n U) i 1 i 2 ...i n-1 j n i n+1 ...i N = i n c i 1 i 2 ...i n-1 i n i n+1 ...i N u j n i n ( 42
)
In this work we have employed the PARAFAC (Parallel Factors) decomposition, also known as CANDECOMP or Canonical Decomposition, provided by the MATLAB tensor toolbox [START_REF] Bader | Efficient MATLAB computations with sparse and factored tensors[END_REF]. This technique provides a diagonal tensor R, suitable for the purpose of this work. The r most relevant entries of this tensor are employed in the decomposition, in order to accomplish with a user-predefined tolerance in the representation of the function . To analyze the behavior of the technique, we have extended the previous problem to threedimensional settings u = 0 in (-1, 1) 3 (43)

subjected to the boundary condition u = x 2 + y 2 + z 2 on *

The higher-order SVD decomposition of the function up to the tolerance 10 -6 provides the functions depicted in Figure 11(a)-(c). Only four functions along each spatial direction were necessary.

With these functions thus calculated, we let the method run to obtain the solution to the problem as depicted along each direction in Figure 12. The global, three-dimensional, solution to the problem is depicted in Figure 13. In this case, 18 functions were necessary along each spatial direction, that is a sum decompositions consisting of 18 sums, to reach the before-mentioned accuracy. This figure also depicts the trace of the solution at the plane x = 0.

A numerical model involving 10 9 degrees of freedom

In order to show the potential of the proposed method, we have analyzed the problem The analytical solution to this problem is given by the function u = (1x 2 )(1y 4 )(1z 6 ). The solution has been approximated by using a separated representation where each one-dimensional function along each coordinate was approximated by using standard piecewise linear FE shape functions on a grid consisting of 1001 nodes on each direction.

u = 2(1-y 4 )(1-z 6 )+12(1-x 2 )y 2 (1-z 6 )+30(1-x 2 )(1-y 4 )z 4 in (-1, 1)
The different computed one-dimensional functions are depicted in Figure 14(a)-(c). As expected an accurate enough solution was obtained after only one iteration, as was the case 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], because the particular form of the exact solution that consists of a single functional product.

The solution obtained is represented in Figure 15, where derefinement has been applied to the results by employing a coarser mesh (for obvious reasons, it is difficult to depict the results on a mesh of 10 9 nodes). The error obtained is depicted in Figure 16, where a slice on an arbitrary section of the domain has been performed for clarity. In this example we do not pursue accuracy, but testing the possibility of employing a huge number of degrees of freedom. What is noticeable about this solution, which approaches the limits of traditional FEs, is that, despite the extremely high number of degrees of freedom involved, the Matlab code ran in around 10 s on a PC with four processors (only one was employed, no parallel calculation was allowed, even if parallelization is an easy matter in the context of separated representations) with 1.96 Gb of RAM memory each, running under Linux.

Models defined in domains other than hypercubes

To show that the proposed method is applicable to problems defined over irregularly shaped domains, other than hypercubes, we consider the following problem defined over a cylindrical domain whose geometry is depicted in Figure 17. The problem is defined as u = 1.0, with u = 0 on * .

Since the analytical solution is not known, we have compared the results with those from a finite element simulation over a fine enough mesh, using standard linear triangular elements over the actual geometry of the domain. By applying the technique explained in Section 3.1, we computed the characteristic function depicted in Figure 18.

In order to apply the separated representation strategy, we first embed the geometry in a squared in which each axis contains 60×60 nodes. The solution outside the domain of interest depicted in Figure 17 is obviously discarded. By comparing the separated representation solution to that obtained by a dense enough FE mesh, and by using the norm e 2 = (1/n) n i e 2 i (with n the number of points in which the error is evaluated and e i the pointwise error at that points), we have obtained after the first iteration an error of 2.97×10 -5 , with no significant improvement by adding other sums of the decomposition. The computed solution is depicted in Figure 19. If more accurate solutions are desired, higher nodal densities along each direction should be employed. It is important to note that the separated representation solution cannot produce solution more accurate than the ones related to a finite element FE discretization on the mesh resulting from the one-dimensional nodal distributions tensor product (Figure 20).

CONCLUSIONS

Problems defined in high-dimensional spaces are very common in Engineering and Applied Sciences. The common feature of all these models relies in the exponentially growing size of the mesh with the number of dimensions of the problem. In previous works by the authors, a technique based on the use of a separated representation consisting of a finite sums decomposition has proven to be very useful for some of these problems.

However, the theoretical establishment of these techniques still deserves further developments. In this paper we have dealt with two particular features of the method, never until now, in our knowledge, addressed. The first one is the imposition of non-homogeneous essential boundary conditions and the second one is related to the simulation in domains other than hypercubes.

The proposed technique is based upon the construction of a change of variable for the first case. This change of variable must be compatible, of course, with the separated representation strategy needing the use of multidimensional SVD techniques for that purpose. The second issue related to the treatment of complex geometries has been accomplished through the method of the characteristic function, a very classical method, recently recovered within the framework of meshless methods. Using both strategies we have solved some numerical models proving the usefulness of the proposed techniques and, in turn, enlarge its field of application. We believe that the proposed techniques have an enormous potential in a wide range of applications, and this constitutes our current effort of the research.

The extension of these strategies to some physical models defined in high dimensions and the consideration of general boundary conditions including general natural conditions are some works in progress.

Figure 1 .

 1 Figure 1. Definition of the boundary of a given domain by a set of six inequalities.

Figure 2 .

 2 Figure 2. Essential boundary of a given domain (a) and associated R-function (b).

  (a)) by black lines and the resulting R-function is depicted (Figure 2(b)).

Figure 3 .

 3 Figure 3. Function obtained.

Figure 4 .

 4 Figure 4. Function for the second problem, Equation (39).

Figure 5 .Figure 6 .

 56 Figure 5. Singular value decomposition of up to the precision 10 -5.

Figure 7 .Figure 8 .

 78 Figure 7. Solution obtained by standard FEM, assumed as the reference solution of the problem.

Figure 9 .

 9 Figure 9. Solution obtained by the proposed technique. Mesh consisting of 80 nodes in each direction.

Figure 10 .

 10 Figure 10. Solution obtained by the proposed technique for the homogeneous problem in z variable. Mesh of 80 nodes.

Figure 11 .

 11 Figure 11. Functions obtained along x (a), y (b) and z (c) spatial directions for the boundary conditions of the problem stated in Equation (43).

Figure 12 .

 12 Figure 12. Functions obtained along x (a), y (b) and z (c) spatial directions for the problem stated in Equation (43).

Figure 13 .

 13 Figure 13. Solution obtained by the proposed technique for the problem with non-homogeneous boundary conditions in three dimensions. Mesh of 40 elements along each direction.

Figure 14 .

 14 Figure 14. Functions obtained along x (a), y (b) and z (c) spatial directions for the problem stated in Equation (44).

Figure 15 .

 15 Figure 15. Solution of the problem stated in Equation (14).

Figure 16 .

 16 Figure 16. Slice on the error map of the problem stated in Equation (14).

Figure 17 .Figure 18 .

 1718 Figure 17. Geometry of the problem in Section 4.5.

Figure 19 .Figure 20 .

 1920 Figure 19. Solution for the problem in Section 4.5.

NUMERICAL EXAMPLES

In this section we describe some simple numerical experiments for illustrating the potential of the discretization technique based on the use of separated representations.

A first example

Let us consider the problem u = 0 in (0, 1) 3 (36) with u = x on * ≡

In this case, essential boundary conditions can be expressed in the form given by Equation ( 24) as The weak form of the problem looks:

Solving the problem in variable z gives the trivial solution z = 0 and thus we have

which coincides with the analytical solution of the problem, as expected.

A second, more elaborated, problem

Let us consider now the problem

subjected to the boundary condition u = x 2 + y 2 on * that in this case do not verify the PDE (39). In this case, with the help of the approximation proposed in Equation (25), we obtain the following function verifying essential boundary conditions:

=

-2+ x 4 + y 4 -2+ x 2 + y 2 (40)


Amine Ammar

Francisco Chinesta

Elías Cueto

Manuel Doblaré

Proper generalized decomposition of time-multiscale models

Keywords: time multiscale, proper generalized decomposition, separated representations, transient models

Models encountered in computational mechanics could involve many time scales. When these time scales cannot be separated, one must solve the evolution model in the entire time interval by using the finest time step that the model implies. In some cases, the solution procedure becomes cumbersome because of the extremely large number of time steps needed for integrating the evolution model in the whole time interval. In this paper, we considered an alternative approach that lies in separating the time axis (one-dimensional in nature) in a multidimensional time space. Then, for circumventing the resulting curse of dimensionality, the proper generalized decomposition was applied allowing a fast solution with significant computing time savings with respect to a standard incremental integration.

INTRODUCTION

Many problems in science and engineering remain intractable despite the impressive progresses in computer science and the computational resources available today, because their numerical complexity is simply unimaginable. Among the models that remain intractable today, we can distinguish two main families:

The first family of models consists of standard models usually encountered in computational mechanics, defined in large and complex three-dimensional (3D) geometries, involving many multiphysics couplings, many scales (in space and time), strong nonlinearities, and whose transient simulation needs extremely small time steps. These models are usually encountered in the mechanics of structures, but they are also present in many other fields. To illustrate this scenario, one could imagine the simple reaction-diffusion model that describes the degradation of plastic materials, where the characteristic time of the chemical reaction involved in the material degradation is of some microseconds, and the one related to the diffusion of chemical substances (that also represents the material degradation characteristic time itself) is of the order of years. In these cases, standard incremental techniques must be replaced by other more efficient techniques. The other family of challenging models concerns those models defined in highly dimensional spaces. This kind of models appears naturally, for example, in the modeling of the structure and properties of materials at the finest scales. Thus, models in quantum chemistry (e.g., the Schroedinger or Dirac equations), the kinetic theory description of materials (e.g., the Fokker-Planck equation), models in financial mathematics, genetics (the chemical master equation) . . . are defined in spaces involving hundreds, thousands, or millions of dimensions. These models exhibit the redoubtable curse of dimensionality when the usual mesh-based discretization techniques are applied. The curse of dimensionality can be easily illustrated. Imagine a model defined in a space involving D dimensions. If one proceeds to the solution of such model 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 D . With M 10 3 (a very coarse description in practice) and D 30 (a quite simple system in practice), the numerical complexity results 10 90 . It is important to recall that 10 80 represents the presumed number of elementary particles in the whole universe. Moreover, as we illustrate later, many usual models can be expressed by introducing many other extra coordinates that will be called configurational or conformational coordinates. For example, model parameters (e.g., the material conductivity involved in a thermal model, the applied loads, the initial or boundary conditions, parameters describing the geometry . . . ) could be considered as extra coordinates. Thus, by solving the resulting multidimensional model, we could have access to the model solution at any point of the space domain, at any time, and for any value of the thermal conductivity, source term, parameter describing the geometry . . . The price to be paid is the solution of a problem that, until now, no numerical technique was able to perform.

Recently, a novel efficient technique able to circumvent the challenging issues just described was proposed. It was coined as proper generalized decomposition (PGD). It is based on a separated representation of the unknown field. Thus, for circumventing the numerical issues associated with the first family of models (previously described), one could consider the transient solution under the separated form:

u.x, t/ N X iD1 X i .x/ T i .t / (1)
representation, which was originally proposed in the 1980s by Pierre Ladeveze as one of the main ingredients of the LATIN (nonlinear and non-incremental solver), and which was called "radial approximation" [START_REF] Ladeveze | Nonlinear Computational Structural Mechanics[END_REF][START_REF] Ladeveze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF][START_REF] Nouy | Multiscale computational strategy with time and space homogenization: a radial-type approximation technique for solving microproblems[END_REF][START_REF] Ladeveze | The LATIN multiscale computational method and the proper generalized decomposition[END_REF][START_REF] Ammar | Non-incremental strategies based on separated representations: applications in computational rheology[END_REF]. If we look at the performance of such separated representation-based discretization techniques, the verdict is for many models implacable. If one considers a standard transient model defined in a 3D physical space, and if one considers P time steps, usual incremental strategies must solve P (in general nonlinear) 3D problems (do not forget that P can be millions). However, if space-time PGD (radial approximation) is considered, we should solve around N m 3D problems for computing the space functions and N m one-dimensional (1D) problems for computing the time functions (m being the number of iterations needed for computing each term of the finite sum because of the nonlinear nature of the definition of the PGD). As m 10 and N 10 for many models, the computing time savings can reach many orders of magnitude.

If we come back to the second family of models, an appealing choice consists of expressing the unknown field as a finite sum of functional products, that is, expressing a generic multidimensional function as

u.x 1 , , x N / N X iD1 X 1 i .x 1 / : : : X D i .x D / (2)
where x i could represent any coordinate, scalar or vector, involving the physical space, the time, or any other conformation coordinate. It could also represent random variables for the modeling of parametric uncertainties. This separated representation was proposed by A. Ammar and F. Chinesta some years ago in the context of multidimensional models encountered in the description of complex fluids within the kinetic theory framework [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 representations[END_REF], as well as by A. Nouy in the context of stochastic models [START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF]. The reader can refer to [START_REF] Chinesta | Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF] and the references therein for some more recent contributions on this topic.

Within the PGD framework for solving a problem defined in a space of dimension D, if M nodes are used to discretize each coordinate, the total number of unknowns involved in the solution is M D instead of the M D degrees of freedom involved in mesh based discretizations. We must recall that these functions are not known "a priori" but are computed by introducing the approximate separated representation into the model weak form and then by solving the resulting nonlinear problem associated with the construction of the separated representation. As it can be noticed in the expression of the separated representation, the complexity scales linearly with the dimension of the space in which the model is defined, instead of the exponential growing characteristic of mesh-based discretization strategies. In general, for many models, the number of terms N in the finite sum is quite reduced (a few tens).

We would like to recall that, using these approaches, we have reduced the computing time related to the solution of problems belonging to the first family of models (defined earlier) in several orders of magnitude, and on the other hand, we solved successfully highly multidimensional models because they were considered suffering of the irremediable curse of dimensionality [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF][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].

Because of the novelty and the youth of the PGD method, many aspects have not yet been addressed. The aim of this work is to push back the limits of these methods. Definitively, these developments could lead to a real change of paradigm in computational mechanics. Imagine the possibility of solving efficiently and with a controlled accuracy any parametric (including parameterized geometries) multiscale and multiphysics model. Inverse identification and optimization would be a simple post-treatment.

In some of our recent works, we focused on multiphysics models involving different characteristic times and behaviors (local and global). This problem was addressed in particular in [START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF]. In the present work, we focus in a slightly different issue, the one related to models involving nonseparable time scales that should be solved by assuming the finest time step concerned by the model solution.

Time multiscale problems involve similar difficulties than those models that are multiscale in space but add some extra issues that must be solved with caution. This is partly caused by the inherent sequential nature of time, which is not present in standard multiscale models. Many works have been devoted in the last years to the topic of establishing efficient time multiscale numerical methods. These approaches are also very different in nature.

For instance, in [START_REF] Yu | Temporal homogenization of viscoelastic and viscoplastic solids subjected to locally periodic loading[END_REF], a method is developed that applies the same principles of spatial homogenization to time multiscale problems. To that end, a neat separation of time scales is mandatory, as in traditional homogenization approaches.

Particularly noteworthy is the topic of establishing parallel approaches in time, needed for computational efficiency, or simply when different physics with different intrinsic time scales are considered. Here, again the sequential nature of time plays an important role and poses important difficulties. [START_REF] Frantziskonis | Time-parallel multiscale/multiphysics framework[END_REF] and [START_REF] Baffico | Parallel-in-time molecular-dynamics simulations[END_REF] are two particular instances of the application of the so-called time-parallel approaches to time-multiscale problems. In [START_REF] Bottasso | Multiscale temporal integration[END_REF], however, the principles of the variational multiscale method are applied to time multiscale integration. In [START_REF] Neron | A computational strategy for poroelastic problems with a time interface between coupled physics[END_REF][START_REF] Neron | A computational for thermo-poroelastic structures with a time-space interface coupling[END_REF], on the contrary, a method based on the LATIN approach is developed that couples different physical problems with different intrinsic time scales in a space-time framework that can be considered as the origin of the approach presented herein.

We consider, however, an alternative approach that lies in separating the time axis (onedimensional in nature) in a multidimensional time space. Then, for circumventing the resulting curse of dimensionality, the PGD is applied allowing for a fast solution with significant computing time savings with respect to a standard incremental integration.

In the next section, we revisit the main ideas of the PGD whose implementation is illustrated through the solution of an academic problem. Then, a tensor formulation of the PGD is presented in Section 3. Section 4 concerns the introduction of constraints in the solution process from the use of penalty or Lagrange multiplier techniques. Section 5 makes use of this procedure for decomposing the time axis into a two-dimensional (2D) time domain procedure that is applied in a quite simple transient ODE for illustrating the capabilities of the proposed technique. In that section, we consider also the solution of a transient parabolic PDE in which the time coordinate is transformed into a 2D time space. In the last case, the resulting model includes a space coordinate and two time coordinates. Finally, Section 6 presents an alternative formulation for enforcing the solution continuity in time without using either penalty or Lagrange multipliers. The continuity is enforced a posteriori.

THE PROPER GENERALIZED DECOMPOSITION REVISITED

In what follows, the construction of the PGD is illustrated by considering the parametric heat transfer equation:

@u @t ku f D 0 (3)
where .x, t , k/ 2 I = and, for the sake of simplicity, the source term is assumed constant, that is, f D cst. Because the conductivity is considered unknown, it is assumed as a new coordinate defined in the interval =. Thus, instead of solving the thermal model for different values of the conductivity parameter, we prefer introducing it as a new coordinate. The price to pay is the increase of the model dimensionality; however, as the complexity of PGD scales linearly with the space dimension, the consideration of the conductivity as a new coordinate remains compatible with fast and cheap solutions.

The solution of Equation ( 3) is sought under the form:

u .x, t , k/ iDN X iD1 X i .x/ T i .t / K i .k/ (4)
In what follows, we are assuming that the approximation at iteration n is already achieved:

u n .x, t , k/ D iDn X iD1 X i .x/ T i .t / K i .k/ (5)
and that, at current iteration, we look for the next functional product X nC1 .x/ T nC1 .t / K nC1 .k/ that, for alleviating the notation, will be denoted by R . We consider the global weak form of Equation ( 3):

Z I = u  @u @t ku f à d x dt d k D 0 (6)
where the trial and test functions write, respectively:

u .x, t , k/ D iDn X iD1 X i .x/ T i .t / K i .k/ C R .x/ S .t / W .k/ (7)
and

u .x, t , k/ D R .x/ S .t / W .k/ (8)
Introducing Equations ( 7) and ( 8) into Equation ( 6), it results in

Z I = R S W Â R @S @t W k R S W Ã d x dt d k D D Z I = R S W R n d x dt d k (9)
where R n defines the residual at iteration n that reads:

R n D iDn X iD1 X i @T i @t K i iDn X iD1 k X i T i K i . f (10)
Now, knowing all the functions involving the time and the parametric coordinate, we can integrate Equation ( 9) in its respective domains I =. Integrating in I = and taking into account the following notations 2 6 6 6 6 6 6 6 6 6 6 4

w 1 D R = W 2 d k s 1 D R I S 2 dt r 1 D R R 2 d x w 2 D R = kW 2 d k s 2 D R I S dS dt dt r 2 D R R R d x w 3 D R = W dk s 3 D R I S dt r 3 D R R dx w i 4 D R = W K i d k s i 4 D R I S d T i dt dt r i 4 D R R X i d x w i 5 D R = kW K i d k s i 5 D R I S T i dt r i 5 D R R X i d x 3 7 7 7 7 7 7 7 7 7 7 5 (11)
Equation ( 9) reduces to:

Z R .w 1 s 2 R w 2 s 1 R/ d x D D Z R iDn X iD1 w i 4 s i 4 X i iDn X iD1 w i 5 s i 5 X i w 3 s 3 f ! d x (12)
Equation ( 12) defines an elliptic steady state boundary value problem that can be solved by using any discretization technique operating on the model weak form (finite elements, finite volumes, . . . ). Another possibility consists in coming back to the strong form of Equation ( 12):

w 1 s 2 R w 2 s 1 R D D iDn X iD1 w i 4 s i 4 X i iDn X iD1 w i 5 s i 5 X i w 3 s 3 f ! (13)
that could be solved by using any collocation technique (finite differences, SPH . . . ).

Computing S .t / from R .x/ and W .k/:

In the present case, the test function writes:

u .x, t , k/ D S .t / R .x/ W .k/ (14)
Now, the weak form reads Z

I = S R W Â R @S @t W k R S W Ã d x dt d k D D Z I = S R W R n d x dt d k (15)
that integrating in the space = and taking into account the notation (11) results in:

Z I S Â w 1 r 1 dS dt w 2 r 2 S Ã dt D D Z I S iDn X iD1 w i 4 r i 5 d T i dt iDn X iD1 w i 5 r i 4 T i w 3 r 3 f ! dt (16)
Equation ( 16) represents the weak form of the ODE defining the time evolution of the field S that can be solved by using any stabilized discretization technique (SU, Discontinuous Galerkin, . . . ).

The strong form of Equation ( 16) reads:

w 1 r 1 dS dt w 2 r 2 S D D iDn X iD1 w i 4 r i 5 d T i dt iDn X iD1 w i 5 r i 4 T i w 3 r 3 f ! ( 17
)
than can be solved by using backward finite differences, or higher order Runge-Kutta schemes, among many other possibilities.

Computing W .k/ from R .x/ and S .t /:

In the present case, the test function writes:

u .x, t , k/ D W .k/ R .x/ S .t/ (18)
Now, the weak form reads

Z I = W R S Â R @S @t W k R S W Ã d x dt d k D D Z I = W R S R n d x dt d k (19)
that integrating in the space I and taking into account the notation [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] results in:

Z = W .r 1 s 2 W r 2 s 1 k W / d k D D Z = W iDn X iD1 r i 5 s i 4 K i iDn X iD1 r i 4 s i 5 k K i r 3 s 3 f ! d k (20)
Equation ( 20) does not involve any differential operator. The strong form of Equation ( 20) reads:

.r 1 s 2 r 2 s 1 k/ W D iDn X iD1 r i 5 s i 4 r i 4 s i 5 k K i r 3 s 3 f ! (21)
that represents an algebraic equation. Thus, the introduction of parameters such as additional model coordinates has no noticeable effect in the computational cost, because the original equation does not contain derivatives with respect to those parameters. Finally, note that other minimization strategies have been proposed, leading to more robust and faster convergence for building up the PGD [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF]. The issues related to nonlinear behaviors, as for example, a thermal conductivity depending on the temperature field, was deeply considered in [START_REF] Ammar | Non-incremental strategies based on separated representations: applications in computational rheology[END_REF] and [START_REF] Pruliere | On the deterministic solution of multidimensional parametric models using the proper generalized decomposition[END_REF]. On the other hand, the issue related to the enforcement of nonhomogeneous boundary conditions was treated in [START_REF] Gonzalez | Recent advances on the use of separated representations[END_REF].

Remark 1

The construction of each term in the sum (4) needs a certain number of iterations because of the nonlinearity of the problem related with the approximation (5). We denote by m i the number of iterations that were needed for computing the i-sum in Equation (4). Let m D P iDN iD1 m i be the total number of iterations involved in the construction of the separated approximation (4). It is easy to note that the solution procedure needs the solution of m 3D problems related to the construction of the space functions X i .x/, i D 1,

, N , m 1D ODEs related to the construction of functions T i .t / and m linear systems related to the definition of functions K i .k/. In general, m i rarely exceeds ten. On the other hand, the number of terms in the sum N needed to approximate the solution of a given problem depends on the solution regularity itself, but all the experiments carried out until now reveal that this number ranges from a few tens to a few hundreds. Thus, we can conclude that the complexity of the solution procedure is of some hundreds of 3D solutions (the cost related to the 1D problems being negligible with respect to the one related to the 3D problems). Now, if we assume a classical approach, one should solve a 3D problem at each time step and for each value of the parameter k. In usual applications, the complexity reaches millions of 3D solutions. In [START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF], we proved that the CPU time savings, by applying the PGD, can be of several orders of magnitude.

Remark 2

We noticed that, depending on the initial choice of functions R.x/, S.t/, and W .k/, the nonlinear iteration algorithm converges to different functions, but when they were normalized by dividing all of them by their respective norms, the resulting functions where the same independently of the initial choice.

PROPER GENERALIZED DECOMPOSITION TENSOR FORM

The procedure described in Section 2 can be generalized by using a tensor notation. We assume that the discrete problem that we write formally as U T AU D U T B, can be expressed in a separated form:

A D n A X iD1 A i 1 ˝Ai 2 ˝ ˝Ai D B D n B X iD1 B i 1 ˝Bi 2 ˝ ˝Bi D U N X iD1 u i 1 ˝ui 2 ˝ ˝ui D ( 22
)
where A i , B i , and u i involve only the coordinate x i .

For an efficient application of the PGD, the differential operators must be separated. A separated representation can be used even if the operators are not separable, but in that case, the integration of the week form becomes cumbersome.

The separated representation of A and B comes directly from the differential operators involved in the PDE weak form.

At iteration n, vectors u i j , 8i Ä n, and 8j Ä D are assumed known. Now, we are looking for an enrichment:

U D n X iD1 u i 1 ˝ ˝ui D C R 1 ˝ ˝RD (23)
where R i , i D 1, , D, are the unknown enrichment vectors. We assume the following form of the test field:

U D R 1 ˝R2 ˝ ˝RD C : : : C R 1 ˝ ˝RD 1 ˝R D (24)
Introducing the enriched approximation into the weak form, the following discrete form results:

n A X iD1 n X j D1 .R 1 / T A i 1 u j 1 .R D / T A i D u j D C C C n A X iD1 n X j D1 .R 1 / T A i 1 u j 1 .R D / T A i D u j D C C n A X iD1 .R 1 / T A i 1 R 1 .R D / T A i D R D C C C n A X iD1 .R 1 / T A i 1 R 1 .R D / T A i D R D D D n B X iD1 .R 1 / T B i 1 .R D / T B i D C C .R 1 / T B i 1 .R D / T B i D (25)
For alleviating the notation, we define:

n C X iD1 C i 1 ˝ ˝C i D D n B X iD1 B i 1 ˝ ˝Bi D n A X iD1 n X j D1 A i 1 u j 1 ˝ ˝Ai D u j D ( 26
)
where n C D n B C n A n. This sum only contains known fields. Thus, Equation (25) can be written as:

n A X iD1 .R 1 / T A i 1 R 1 .R D / T A i D R D C C C n A X iD1 .R 1 / T A i 1 R 1 .R D / T A i D R D D D n C X iD1 .R 1 / T C i 1 .R D / T C i D C C .R 1 / T C i 1 .R D / T C i D ( 27
)
This problem is strongly nonlinear. To solve it, a method of alternated directions can be applied. The idea is, starting with the trial vectors R .0/ i , i D 1, , D or assuming known these vectors at iteration r 1, R .r 1/ i , i D 1, , D, update them using an appropriate strategy. The simplest alternatives consist of:

Update vectors R .r/ i , 8i , from R .r 1/ 1 , , R .r 1/ i 1 , R .r 1/ iC1 , , R .r 1/ D . Update vectors R .r/ i , 8i , from R .r/ 1 , , R .r/ i 1 , R .r 1/ iC1 , , R .r 1/ D .
The last strategy converges faster, but the advantage of the first one is the possibility of updating each vector simultaneously making use of a parallel computing platform. The fixed point of this iteration algorithm allows defining the enrichment vectors

u nC1 i D R i , i D 1,
, D. When we look for vector R k , assuming known all the others R i , i ¤ k, the test field reduces to:

U T D R 1 ˝ ˝Rk 1 ˝R k ˝RkC1 ˝RD (28)
The resulting discrete weak form writes:

n A X iD1 R T 1 A i 1 R 1 R T k A i k R k R T D A i D R D Á D D n C X iD1 R T 1 C i 1 R T k C i k R T D C i D ( 29
)
Making use of the arbitrariness of R K , the following linear system can be easily obtained:

0 @ n A X iD1 0 @ D Y j D1,j ¤k R T j A i j R j 1 A A i k 1 A R k D n C X iD1 0 @ D Y j D1,j ¤k R T j C i j 1 A C i k (30)
which can be easily solved.

INTRODUCING SOLUTION CONSTRAINTS

In numerical homogenization techniques, a classical procedure consists of considering different time scales as different dimensions of the problem. In principle, this seems to be specially well suited for a PGD treatment, especially if the number of scales is high. But the question is, now if the PGD can effectively treat the same problem, is there no neat separation of scales? Converting the time axis in a multidimensional space forces us to enforce the continuity on the resulting domain boundary. This section develops the PGD formulation in the presence of constraints and their tensor form, which are particularly convenient for code implementation.

Using the previous notation, we assume that the problem to be solved writes:

A U D B (31) with A D n A X j D1 A j 1 ˝Aj 2 ˝: : : ˝Aj D (32) B D n B X j D1 B j 1 ˝Bj 2 ˝: : : ˝Bj D ( 33
)
being the solution tensor form:

U n u X j D1 u j 1 ˝uj 2 ˝: : : ˝uj D ( 34
)
The sizes of the different vectors are r 1 , r 2 , .., r D , respectively, whereas the size of matrix

A j i is r i r i , i D 1,
, D. We assume in this section that the solution is searched under the constraints given by:

H U D J (35) with H D n H X j D1 H j 1 ˝H j 2 ˝: : : ˝H j D (36) J D n J X j D1 J j 1 ˝J j 2 ˝: : : ˝J j D ( 37
)
whose particular form will be illustrated later. In the last expression, the sizes of vectors J j i are assumed to be identical, s, being the sizes of matrices H j i are s r i with i D 1, , D.

Penalty formulation

We introduce the functional to be minimized as

E D kA U Bk 2 C ˛kH U J k 2 (38)
with ˛being the penalty parameter.

The previous equation can be rewritten as:

E D .A U B/ T .A U B/ C ˛.H U J / T .H U J / (39)
whose minimization with respect to the unknown field leads to:

dE d U D 0 ) A T A U A T B C ˛ H T H U H T J D 0 (40)
Finally, the system to be solved writes:

A T A C ˛HT H U D A T B C ˛HT J . ( 41
)

Remark 3

This strategy deserves the following comments:

The number of operators involved by

A T A C ˛HT H results in n 2 A C n 2 H . All the terms have dimensions .r 1 , r 1 / ˝: : : ˝.r D , r D /.
The best choice of the penalty parameter is far to be simple. The resulting system is symmetric, making possible the use of solution strategies like the one described in Section 2.

Lagrange multipliers

By introducing the Lagrange multiplier the system to be solved writes:

 A H T H 0 à " ƒ' … K  U à D  B J à " ƒ' … F ( 42
)
where the size of K is .r 1 C s, r 1 C s/ ˝ ˝.r D C s, r D C s/ that can be written as:

K D n A X j D1  A j 1 0 .r 1 ,s/ 0 .s,r 1 / 0 .s,s/ à ˝ ˝Â A j D 0 .r D ,s/ 0 .s,r D / 0 .s,s/ à C n H X j D1  0 .r 1 ,r 1 / H jT 1 0 .s,r 1 / 0 .s,s/ à ˝ ˝Â 0 .r D ,r D / H jT D 0 .s,r D / 0 .s,s/ à C n H X j D1  0 .r 1 ,r 1 / 0 .r D ,s/ H j 1 0 .s,s/ à ˝ ˝Â 0 .r D ,r D / 0 .r D ,s/ H j D 0 .s,s/ à (43)
expression that involves n A C 2n H terms.

Concerning the right hand member, we have

F D n B X j D1  B j 1 0 .s,1/ à ˝ ˝Â B j D 0 .s,1/ à C n J X j D1  0 .r 1 ,1/ J j 1 à ˝ ˝Â 0 .r D ,1/ J j D à . ( 44
)
The searched solution writes in the present case:

U D Â U Ã D n F X j D1 F j 1 ˝F j 2 ˝ ˝F j D (45)
which involves a tensor product of vectors of size .r 1 C s, 1/ ˝ ˝.r D C s, 1/ from which we can easily extract the problem solution U from

U D n u X j D1 u j 1 ˝ ˝uj D (46)
where vector u i contains the first r i entries of vector F i . Representation (46) defines a non-optimal decomposition that could be reduced by applying a multidimensional singular value decomposition [START_REF] Chinesta | Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models[END_REF].

Remark 4

This strategy deserves the following comments:

The main drawback of this approach lies in the fact that the resulting system is nonsymmetric.

It could be symmetrized by minimizing the L 2 norm of the residual:

E D KU F 2 (47)
whose minimization leads to:

dE d U D 0 (48)
that is equivalent to:

K T KU D K T F (49)
The number of terms involved in

K T K is .n A C 2n H / 2 .
The main advantage of this approach is the absence of any adjustable parameter.

Illustrating the prescription of constraints through a simple example

In this section, we are defining the constraints tensor form related to a quite simple example that will be used later for separating the time axis.

The generic tensor form related to the solution constraints writes, using the notation just introduced, as

H U D J (50)
where vectors J j i have the same size s and matrices H j i are of size s r i with i D 1, , D. We assume the situation depicted in Figure 1 that represents a 2D model defined in a square domain whose axes x 1 and x 2 were discretized by using three and f our nodes, respectively, that is, r 1 D 3 and r 2 D 4.

We would like to prescribe the following constraints in the discrete model:

Prescribe the value u d at the node symbolized by a square in Figure 1.

Prescribe the same value of the solution at nodes symbolized in Figure 1 by a triangle. Prescribe the same value of the solution at nodes symbolized in Figure 1 by a star.

These three conditions write: whose matrix form reads as:

® 1 0 0 ˝ 1 0 0 0 ¯U D u d ˝1 (51)
® 1 0 0 ˝ 0 0 0 1 C 0 1 0 ˝ 1 0 0 0 ¯U D 0 ˝0 (52) ® 0 1 0 ˝ 0 0 0 1 C 0 0 1 ˝ 1 0 0 0 ¯U D 0 ˝0 (53)
H D 0 @ 1 0 0 0 0 0 0 0 0 1 A ˝0 @ 1 0 0 0 0 0 0 0 0 0 0 0 1 A C 0 @ 0 0 0 1 0 0 0 0 0 1 A ˝0 @ 0 0 0 0 0 0 0 1 0 0 0 0 1 A C 0 @ 0 0 0 0 1 0 0 0 0 1 A ˝0 @ 0 0 0 0 1 0 0 0 0 0 0 0 1 A C 0 @ 0 0 0 0 0 0 0 1 0 1 A ˝0 @ 0 0 0 0 0 0 0 0 0 0 0 1 1 A C 0 @ 0 0 0 0 0 0 0 0 1 1 A ˝0 @ 0 0 0 0 0 0 0 0 1 0 0 0 1 A (54)
The right-hand member results in:

J D 0 @ u d 0 0 1 A ˝0 @ 1 0 0 1 A (55)

TIME-MULTISCALE PGD FORMULATION WITHOUT SEPARATION OF SCALES

In this section, we are applying the procedure described in the previous section for separating the time axis, firstly, from 1D in nature into a 2D domain, and then to more complex scenarios.

Solving an ordinary differential equation by increasing the temporal dimensionality

For this purpose, we should define the structure of the discrete problem expressed by:

A U D B ( 56
)
related to the transient ODE:

du dt D f .t/ (57)
that we would integrate in the whole time interval using r 1 .r 2 1/ time steps. When the number of time steps is excessive, we could divide the whole time interval OE0, t max into r 1 intervals h 0, t max r 1 i h

t max .r 1 1/ r 1 , t max i .
It is important to note that the end of each one of these intervals corresponds to the beginning of the next one. Thus, we define two dimensions, the first one discrete

x 1 D 0, 1, , .r 1 1/ (58)
and the second one continuous

x 2 D Ä 0, t max r 1 . (59)
Now, the transient equation can be rewritten in the 2D time space as:

@ @x 2 u.x 1 , x 2 / D f .x 1 , x 2 / (60)
The discrete model writes:

A D 0 B B B @ 1 0 0 0 1 0 . . . . . . . . . . . . 0 0 : : : 1 1 C C C A .r 1 ,r 1 / ˝1 t 0 B B B @ 1 1 1 . . . . . . 1 1 1 C C C A .r 2 ,r 2 / (61) with t D t max r 1 .r 2 1/ (62)
and

B D svd .f .x 1 , x 2 // (63)
were the singular value decomposition is applied after discretizing the continuous coordinate x 2 .

Obviously, the separation of f .x 1 , x 2 / is a delicate point. There are different scenarios: (1) proceed to a continuous "singular value decomposition" separation when it is possible as described in [START_REF] Beylkin | Algorithms for numerical analysis in high dimensions[END_REF];

(2) use the coarsest discrete description and then project on the finest mesh; and (3) in some cases, the source term comes from the solution of other equation (in coupled models), and in that case, it has the appropriate separated form.

In what follows, we consider the source term of the transient model given by f .t/ D 2t cos 2 .!t / 2t 2 ! cos.!t / sin.!t / (64)

with ! D 10 and t max D 5. The solution is performed by assuming r 1 D 10 (that constitutes the first, coarse, time coordinate that will be designated by t ) intervals, each one involving r 2 D 41 time steps (that constitute the second time coordinate, the finest one, designated by ). The solution is then searched in the separated form:

u.t , / iDN X iD1 F i 1 .t / F i 2 . / (65)
by enforcing both the verification of the equation as well as the continuity between the end of one interval and the beginning of the next one. This enforcement, either by penalization or by Lagrange multipliers, was deeply addressed in the previous section. Note that, in the present approach, no neat separation of scales is assumed, in deep contrast with the time homogenization approach in [START_REF] Yu | Temporal homogenization of viscoelastic and viscoplastic solids subjected to locally periodic loading[END_REF]. In that work that assumes an additive (instead of multiplicative) enrichment of the macro time scale, the response fields are assumed to be periodic with respect to the fine scale, which is not the case here. Note also that in [START_REF] Neron | A computational strategy for poroelastic problems with a time interface between coupled physics[END_REF] a POD-like approach is employed, in which an averaged fine-scale time evolution is employed, leading to discontinuous time evolutions. This makes it necessary to employ discontinuous Galerkin approaches in time. In the approach here presented, the continuity between time scales is imposed explicitly, as commented before.

Figure 2 depicts the different functions F i 1 .t / and F i 2 . / involved in the solution computed by applying the penalty enforcement of the continuity conditions. Figure 3 shows the computed solution u.t , /. Similar results were obtained when applying the Lagrange multiplier strategy. Finally, both 2D solutions (the one computed by applying the penalty strategy and the one obtained by using Lagrange multipliers) can be pushed down for defining the associated 1D solutions.

When the continuity is enforced by using a penalty technique, the computed solution for different number of terms in the finite sum (65) is depicted in Figure 4. Figure 5 depicts similar results when the continuity is enforced by using lagrange multipliers. We can notice in these figures that by using 10 terms in the sum (65), the separated representation solutions agree in minute to the reference ones computed by using a standard fully incremental integration.

Figure 6 compares the convergence obtained by using both the Lagrange multipliers and the penalty strategies. We can notice that, as expected, the convergence is better when applying Lagrange multipliers. On the other hand Figure 7 compare the computing time needed for constructing both separated representations, and proves that Lagrange multipliers is more expensive from the computational time viewpoint.

We can compare the complexities related to a fully incremental integration and the one related to the construction of separated representation. In the fully incremental integration, we should compute the solution at each one of the r 1 .r 2 1/ time steps that quantify the solution complexity. When a separated representation is built, one should compute N m ODE integrations (m being the iterations involved in the nonlinear solver used for computing each functional product in Equation ( 65)) in the intervals involving r 2 1 time steps. Thus, the complexity results N m .r 2 1/. By comparing both complexities, we can conclude that if N m < r 1 the time axis separation could be an appealing solution for speeding up the integration of transient models. In general, N is of the order of few tens and m is lower than 10. This implies that r 1 should be greater than 100 to induce CPU time savings. However, in this paper, we are more interested in analyzing the possible multidimensional description of the time axis than in proving its numerical advantages.

Solving time-multiscale PDE by increasing the temporal dimensionality

Until now, our attention focused in the solution of time-multiscale ODE. In the present section, we are considering models involving time and space coordinates (PDEs) and exploring the applicability of a decomposition of the time axis in a higher dimensional time space.

For this purpose, we consider the transient PDE equation

@u @t D @ 2 u @x 2 C f .x, t/, .x, t/ 2 .0, / .0, t max (66)
where

f .x, t/ D 2t cos 2 .!t / 2t 2 ! cos.!t / sin.!t / C t 2 cos 2 .!t / sin.x/ (67)
with u.x D 0, t/ D u.x D , t/ D 0, and u.x, t D 0/ D 0.

The exact solution of Equation (66) writes:

u ex .x, t/ D t 2 cos 2 .!t / sin.x/. (68)
For its numerical solution, we make a partition of the whole time interval OE0, t max into r 1 intervals h 0, t max r 1 i h

t max .r 1 1/ r 1 , t max i .
In what follows, we consider t max D 5, r 1 D 10, and ! D 10. Thus, we define two time dimensions, the first one discrete

x 1 D 0, 1, , .r 1 1/ (69)
and the second one continuous

x 2 D Ä 0, t max r 1 . (70)
Moreover, the space dimension x 3 , the present case and without loss of generality, 1D, is discretized by using r 3 nodes. Now, the transient equation can be rewritten in the resulting 3D space as:

@u.x 1 , x 2 , x 3 / @x 2 D @ 2 u.x 1 , x 2 , x 3 / @x 2 3 C f .x 1 , x 2 , x 3 / (71)
For the solution, we consider the Lagrange multiplier strategy for enforcing the continuity with respect to the time coordinate that considers a Lagrange multiplier for each node of the spatial discretization.

Figures 8 and9 depict the time and space functions respectively involved in the separated representation. Figure 10 shows the reconstructed space-time solution post-treated from the 3D separated representation, where the two time coordinates were pushed down for defining a single 1D time axis. This solution is compared with the exact solution, and, as expected, they are in perfect agreement. In order to appreciate the solution features, Figure 11 depicts another view of the computed solution in order to emphasize the presence of fast oscillations.

AN ALTERNATIVE FORMULATION

When the number of dimensions increases, so does the complexity associated with the imposition of time continuity, as discussed in Section 4.3. In this section, an alternative method is developed that has rendered excellent results when applied to ODEs and PDEs.

Application to ODEs

To illustrate the method, we consider firstly the ODE already considered in Section 5.1 @u @t D f .t/ D 2t cos 2 .!t / 2t

Concerning the particular solution, it is sought in the separated form:

u part, compl .t , / Q X j D1 F j t .t / F j . / (75)
Because we are looking for a particular solution, we could enforce it to vanish on the boundary D 0, that is, u part, compl .t , D 0/ D 0. This condition can be obtained by enforcing F j . D 0/ D 0, j D 1, : : : , Q.

The second part of the solution derives from homogeneous equation whose separated form where D t max =r 1 . By using this strategy, the solution of the problem is obtained in only one iteration. The result agrees in minute with the reference one as shown in Figure 12.

u hom .t , / Q 0 X j D1 H j t .
Because of the simplicity of the operators involved in this description, with respect to the Lagrange multiplier strategy discussed previously, the computing time savings in the example just treated is of two orders of magnitude.

Figure 13 compares the evolution of the error with the number of terms considered in the separated representation of the general solution of the homogenous equation Q 0 and the particular solution of the complete equation Q.

For a given accuracy, the number of terms is much lower than that of the one needed when applying the strategies described in Section 5. That is, for the same number of functions in the separated representation of the solution, the one related to the technique described in this section is better than that of the ones computed by using the strategies described in Section 5.

This strategy is generalizable to transient PDEs.

Application to PDEs

For illustrating the application to discretize PDE, we consider the model

@u @t u D f .x, t/ (79)
with, for the sake of simplicity, homogeneous initial and boundary conditions. We consider the model defined in a 2D space domain , x D .x, y/ 2 , with t 2 .0, t max . where the coefficients ˛i j are computed by enforcing the solution time-continuity. Considering a simple nodal collocation for enforcing the time continuity, we can write, for i D 1, , r 1 1:

u part, compl .x k , t i , D / C j DM X j D1 ˛i j u j hom .x k , t i , D / D D j DM X j D1 ˛iC1 j N j .x k / D ˛iC1 k , k D 1, , M (85)
where the fact that N j .x k / D ı jk has been used. Coefficients ˛1 k are computed from the initial condition u.x, t D 0/ D u g .x/:

u g .x k / D u part, compl .x k , t 1 D 0/ C j DM X j D1 ˛1 j u j hom .x k , t 1 , D 0/ D D j DM X j D1 ˛1 j N j .x k / D ˛1 k , k D 1, , M (86)
This technique could be applied in the nonlinear case by performing a linearization prior to proceeding to the solution. However, the main drawback of such an approach lies in the necessity of solving M homogeneous equations. In some applications, M can be extremely high. Obviously, if M is of the same order of magnitude than that of r 1 , the solution procedure has the same complexity than that of a standard incremental solution. For circumventing this difficulty, we propose considering an adaptive reduced basis in order to enforce the initial condition as well as the solution continuity in time.

Thus, we consider two different space meshes: one, fine enough, for approximating the space functions involved in the separated representation F j x,i .x/, and F j x .x/; and the second one, N i .x/, coarser, for enforcing the time-continuity. Because both associated approximation spaces are different by enforcing the time-continuity at the M nodes related to the coarsest mesh, we do not ensure the perfect continuity. The resulting gap could be used to define a criterion for refining the coarsest mesh, the one used to enforce the continuity. However, if one considers standard polynomial approximation bases, as soon as a refinement is performed, the whole solution process must be repeated again. For alleviating the solution process, one could consider a hierarchical basis (very well known in the finite element framework) because in that case, the addition of a new approximation function allows to keep all the work already done. Thus, if we consider the hierarchical basis H i .x/ instead of the classical one N i .x/, with i D 1, : : : , M , if we enrich the discrete approximation space, consisting now in M 0 approximation functions H 0 i .x/, M 0 > M , the M first approximation functions in the refined basis remains unchanged, that is, H 0 i .x/ D H i .x/, i D 1, : : : M . For standard finite element bases, as soon as we change the space dimension, the approximation functions change. Thus, when we proceed to refine a hierarchical approximation basis, the work already done is reused.

When using a hierarchical approximation basis, the solution procedure is slightly modified. If we are considering a discrete hierarchical basis of size M , we should solve the following M homogeneous problems:

@u i hom @t u i hom D 0 u i hom .x, t , D 0/ D H i .x/ μ i D 1, : : : , M (87)
The preceding solutions are sought in the separated form:

u i hom .x, t , / Q 0 X j D1 F j x,i .x/ F j t ,i .t / F j ,i . / (88)
where functions F j x,i are approximated using a fine enough approximation basis. At the first level, we have f ive edge shape functions related to the five internal mesh edges, all of them depicted in Figure 19.

At the second level, we should consider the next five edge shape functions (related to the five internal mesh edges) as well as the six face shape functions, all of them depicted in Figure 20.

Figure 21 depicts the separated representation of the solution u 1 hom of the homogeneous problem:

´@u 1 hom @t u 1 hom D 0 u The separated representation, for the desired precision, involves 15 functional products, that is, 22 when a hierarchical approximation basis consisting of the first two levels (involving 16 approximation functions: f ive first-level edge shape functions, f ive second-level edge shape functions, and finally six second-level face shape functions) was applied.

u 1 hom .x, t
Obviously, the difference between the computed solution and the reference one obtained by using a standard incremental strategy in the whole time interval, decreases as the size of the hierarchical basis increases. Thus, when we use for enforcing the time-continuity only the first-level approximation functions (f ive shape functions) the error is 0.1 (10%) whereas by using the shape functions of the two first levels (16 shape the error reduced to 0.03 (3%). These errors are defined from the L2 norm of the difference between the computed solution and the reference one (the one computed by using a standard incremental strategy) divided by the L2 norm of the reference solution.

CONCLUSION

In this paper, we proved that models involving different non-separable time scales can be reformulated by introducing different time coordinates allowing to express the time dependence of the solution from a multidimensional time approximation. Event if, at present, it is too early for concluding on the computing time and the procedure performances, the possibility of performing such decomposition is, in our opinion, quite interesting, and it could be at the origin of many developments for reducing the computing time of complex multiscale thermomechanical models.

Figure 1 .

 1 Figure 1. 2D problem with constraints.

F2Figure 2 .

 2 Figure 2. Functions involved in the separated representation using a penalty technique for enforcing the continuity constraints: F i 1 .t/ (left) and F i 2 . / (right).

Figure 3 .

 3 Figure 3. Separated representation solution u.t, / computed using a penalty technique for enforcing the continuity constraints.

Figure 4 .

 4 Figure 4. Comparison of the separated solution and the fully incremental one computed for different number of terms N in the separated representation: (top-left) N D 1, (top-right) N D 2, (bottom-left) N D 3, and (bottom-right) N D 10, when a penalty strategy was applied. The green curve represents the approximated solution, whereas the blue one is the reference, exact, solution.

Figure 5 .

 5 Figure 5. Comparison of the separated solution and the fully incremental one computed for different number of terms N in the separated representation: (top-left) N D 1, (top-right) N D 2, (bottom-left) N D 3, and (bottom-right) N D 10, when a Lagrange multiplier was introduced. The green curve represents the approximated solution, whereas the blue one is the reference, exact, solution.

Figure 6 .

 6 Figure 6. Error versus the number of iterations (number of terms in the finite sums decomposition) when using Lagrange multipliers (left) and penalty (right) strategies for enforcing the solution continuity.

Figure 7 .

 7 Figure 7. CPU time versus the number of iterations (number of terms in the finite sums decomposition) when using Lagrange multipliers (left) and penalty (right) strategies for enforcing the solution continuity.

Figure 8 .F3Figure 9 .

 89 Figure 8. Functions of time involved in the separated representation using Lagrange multipliers for enforcing the continuity constraints:F i 1 .x 1 / (left) and F i 2 .x 2 / (right).

Figure 10 .

 10 Figure 10. Space-time reconstructed solution u.x, t/ from the separated representation whose functions were depicted in Figures 8 and 9 (left) and exact solution (right).

Figure 11 .

 11 Figure 11. Another view of the space-time reconstructed solution u.x, t/ depicted in Figure (left).

Figure 12 .

 12 Figure 12. Result for ODE obtained with the technique explained in Section 6.1.

Figure 14 .

 14 Figure 14. Geometry of the domain considered for the solution of transient problem (93).

Figure 15 .Figure 16 .

 1516 Figure 15. Mesh associated to the Szabo and Cernevali hierarchical basis.

Figure 17 .

 17 Figure 17. Edge shape functions related to one of the triangles depicted in Figure 15.

Figure 18 .

 18 Figure 18. Second level edge shape functions and face shape functions related to one of the triangles depicted in Figure 15.

Figure 19 .

 19 Figure 19. Active approximation functions at the first level of the hierarchical approximation related to the mesh depicted in Figure 15.

  1 hom .x, t , D 0/ D H 1 .x/ (94)where H 1 .x/ is the first level edge shape function depicted in Figure19at the top and left.

Figure 20 .

 20 Figure 20. Active approximation functions at the second level of the hierarchical approximation related to the mesh depicted in Figure15:(top) the five second level edge shape functions, and (down) the six first face shape functions.

2 Figure 21 .

 221 Figure 21. Main functions involved in the separated representation of u 1 hom .x, t , /: (top) the first four space functions F i x,1 .x/, i D 1, , 4; (down) the first time functions F i t,1 .t/ (left) and F i ,1 (right), i D 1, , 4.

Figure 22 .

 22 Figure 22. Reconstructed solution u.x, t D 1/ obtained from u.x, t , /.

  D u gen, hom .t / C u part, compl .t / (73) that is, express the solution of Equation (72) as the sum of the general solution of the homogeneous equation plus a particular solution of the complete equation. By performing the time-axis decomposition, it results in u.t , / D u gen, hom .t , / C u part, compl .t , /.

	(72)
	with u.t D 0/ D 0.
	Because of the linearity of Equation (72), an alternative method to find a numerical solution to
	this problem consists in writing
	u.t /

2

! cos.!t / sin.!t /

  t / H j . / (76) that is integrated by enforcing the boundary condition u .t , D 0/ D 1. For this purpose, we prescribe H 1 t .t / D 1 and H 1 . D 0/ D 1 and H j . D 0/ D 0, j D 2, , Q 0 . Thus, the general solution of the homogeneous equation reads: u gen, hom .t i , / D ˛i u hom .t i , / (77) The coefficients ˛i are determined by enforcing the solution time-continuity, by solving sequentially ² ˛1 D u.t D 0/ D 0 u part, compl .t i , / C ˛i u hom .t i , / D ˛iC1 , i D 1, : : : , r 1 1 (78)

  , / significative functions are depicted in Figure 21. We depict, in this figure, the first four space functions F i x,1 .x/, i D 1, , 4, and the first time functions F i t ,1 .t / and F i ,1 , i D 1, , 4. The reconstructed solution u.x, t D 1/ is depicted in Figure

	15 X	F i x,1 .x/ F i t ,1 .t / F i ,1 . /	(95)
	iD1
	whose four most

ACKNOWLEDGEMENT

Contract/grant sponsor: Spanish Ministry of Education and Science; contract/grant number: CICYT-DPI2008-00918

We compute the particular solution of the complete equation in the separated form

with u part, compl .x, t , D 0/ D 0.

The general solution of the homogeneous equation needs solving

for any initial condition.

If we are considering a mesh composed of M nodes for approximating the field u in the space domain , we could solve the following M homogeneous problems

with N i .x/ the elements of a discrete space approximation basis. In the context of finite element techniques, N i .x/ represent the standard shape functions verifying the Kronecker's delta property, that is,

The preceding solutions are sought in the separated form:

With all these solutions u i hom .x, t , /, i D 1, , M , known, the general solution of the homogeneous equation reads:

With all these solutions u i hom .x, t , /, i D 1, , M , known, the general solution of the homogeneous equation reads:

where the coefficients ˛i j are computed by enforcing the solution continuity in time. Considering a simple nodal collocation at M locations x k for enforcing the time-continuity, we can write, for i D 1, , r 1 1:

Coefficients ˛1 k are computed from the initial condition u.x, t D 0/ D u g .x/:

for k D 1, : : : , M . Obviously, by enforcing the solution continuity in time at positions x k , k D 1, , M , we do not ensure the continuity in time everywhere. Thus, the resulting gap can be used as a criterion for enriching the coarse approximation space whose size becomes M 0 > M . However, because its hierarchical nature of the approximation functions, all the already computed solutions u i hom .x, t , /, i D 1, : : : , M , remain unchanged, and then, we only need to compute the M 0 M new homogeneous problems

For illustrating the procedure, we are considering the linear parabolic PDE

with homogeneous initial and boundary conditions. We consider x 2 with depicted in Figure 14, t 2 .0, 1, and ! D 45.

The hierarchical approximation basis consists of the Szabo and Cernevali basis involving vertex, edge, and face shape functions [START_REF] Szabo | Introduction to Finite Element Analysis[END_REF]. The hierarchical basis is associated to the coarse mesh depicted in Figure 15. Figure 16 superpose the fine mesh used for solving the different equations, and the corse mesh used for defining the hierarchical basis employed for enforcing the continuity in time.

The zero level of the hierarchical basis consists on the standard vertex linear finite element shape functions. The first level consists of edge functions. In Figure 17, we depict the ones related to one of the triangles of the mesh shown in Figure 15.

The second level consists of other edge shape function as well as a face function defined on each triangle. Figure 18 depicts those functions for one of the triangles involved in the hierarchical mesh.

Because of the mesh used in our calculations (Figure 15) and the assumed homogeneous boundary conditions, we are not considering the vertex shape functions because all of them are located on the domain boundary @ in which the solution vanishes. By the same reasons, we are not considering all the edge shape functions of any level related to edges located on the domain boundary @ . Thus, we are only retaining the edge shape functions associated with internal edges as well as all the face shape functions.