Alleviating Mesh Constraints: Model Reduction, Parallel Time Integration and High Resolution Homogenization

Industrial processes involving composite materials need for eﬃcient numerical simulations in order to optimize the process parameters. Even if the thermo-mechanical models are nowadays well established, eﬃcient simulations need for further developments. In this work we are addressing some of these issues, in particular the one related to fast solutions combining model reduction and parallel time integration. A separated representation will be also proposed in the context of material homogenization allowing to alleviate the usual mesh constraints.


Introduction
Industrial processes involving composite materials need for efficient numerical simulations in order to optimize the process parameters.Even if the thermo-mechanical models are nowadays well established, efficient simulations need for further developments.
In this work we are considering some of these issues, analyzed from a methodological point of view, without considering its industrial counterpart that requires the coupling of different numerical procedures, some of them proposed and analyzed in this work.
Thermal models involved in the numerical modelling of composite plates welding processes introduce, despite its geometrical simplicity, a certain number of numerical difficulties related to: (i) the very fine mesh required due to the small domain thickness with respect to the other characteristic dimensions as well as to the presence of a heat source moving on the domain surface; (ii) the long simulation times induced by the small thermal conductivity of polymers, the important thermal solicitations as well as to the movement of the heat source; and (iii) the necessity to define a homogenized thermal conductivity which can vary significantly from one point to other in the domain.
The consequence of the first difficulty is the necessity of using very fine meshes everywhere, and therefore, a great number of degrees of freedom that could affect seriously the algorithms efficiency.Moreover, when the transient problem must be solved incrementally, the extremely large number of time steps induces a second difficulty.For alleviating this drawback we will propose in the first part of the paper a parallel time discretization combined with an adaptive reduced modelling.
The other difficulty concerns the material heterogeneity at the microscopic scale.Even if we could in a first approximation neglect the dependence of thermal parameters (specific heat, thermal conductivity, etc.) on the temperature -linear modelling -the macroscopic thermal conductivity will depend on the considered point in the macroscopic part, because it depends strongly on the microstructure details (fibers volume fraction, fiber orientations, etc.).Due to the small size of fibers, very fine models are required in order to analyze the microscopic patterns, which implies the necessity of using high resolution homogenization techniques.As just discussed, the main drawback comes from the extremely large number of degrees of freedom involved in a such microscopic analysis where the level of detail could be in the order of the size of a fiber.In the second part of the present work we will applied a novel numerical technique, recently proposed in [2], that could be an appealing choice for treating this kind of numerical models.
In both topics, the parallel time integration and the homogenization techniques there are some well established approaches that will be discussed later in the text.The main aim of this paper is the coupling of model reduction and time parallel integration in a way never until now considered, at least in our knowledge.The proposed procedure could be considered alternative to the ''parare ´elle'' approach proposed in [14].The second contribution of this work is the application of a separated representation, originally proposed in the context of multidimensional models [2,3], for performing high resolution homogenization.
In all the cases we are only describing the first steps in the development of all these numerical strategies, whose analysis, integration and optimization require further developments that are at present in progress.

Thermal model
From now on, the thermal model is assumed defined in the 3D plate domain X = ]0, L x [ • ]0, L y [ • ]0, L z [ where L z ( L x and L z ( L y , in the time interval [0, t max ], on whose boundary oX different kind of boundary conditions are prescribed, the temperature on o T X and the heat flux on o Q X, with oX = o T X [ o Q X.The components of x 2 X will be denoted by (x, y, z).Without loss of generality, we consider o T X = {x : z = L z } and o Q X the complementary part of the domain boundary.The temperature T(x, t) is then prescribed on o T X to be a given temperature T g (x, t).In o Q X, and for the sake of simplicity in the description of the proposed numerical approaches, we assume a null heat flux across the boundary, that is (K(T, x)Grad T) AE n = 0, where K(T, x) denotes the homogenized thermal conductivity tensor and n the unit outwards vector defined on the domain boundary.Moreover, an initial condition must be imposed at the initial time: T(x, t = 0) = T 0 (x).
The usual variational formulation of the conduction heat equation is expressed by: Find where H 1 is the usual Sobolev functional space, V = {v 2 H 1 : v(x 2 o T X, t) = 0}, and the specific heat and the homogenized thermal conductivity tensor depend on space (in non-homogeneous media) and on the temperature.The source term f(x, T, t) accounts for possible couplings or change of phase phenomena.
From now on, as we consider null thermal fluxes on o Q X the boundary integral in Eq. ( 1) vanishes.Moreover, the thermal properties in the domain X will be considered uniform, implying C = C(T) and K = K(T), which results in To solve Eq. ( 2) an appealing technique lies in the use of the finite element method (FEM) which operates on a mesh of the domain X ¼ [ e¼N e e¼1 X e .Due to the small domain thickness and the evolution of the prescribed temperature on the upper-boundary of the plate, the mesh size must be very fine, with a consequent impact on the number of nodes, and then on the number of degrees of freedom.
Introducing appropriate continuous field approximations in Eq. ( 2), we obtain assuming a semi-implicit time integration: where vector T contains the nodal temperatures (the associated superscripts refer the considered time step) and F accounts for the prescribed nodal temperatures and the thermal sources in X. M and D represent specific heat and thermal conduction effects, respectively.This equation can be rewritten in a more compact form as The main issue related to Eq. ( 4) is the large size of the discrete problem to be solved at each time step.To alleviate this drawback we will consider a model reduction allowing significant CPU time savings.For this purpose we start introducing the proper orthogonal decomposition (POD) also knows as Karhunen-Loe `ve decomposition.

The proper orthogonal decomposition
We assume that the evolution of a certain field T(x, t) is known.In practical applications, this field is defined at the nodes of a spatial mesh x i and for some times T m defines the vector containing the nodal degrees of freedom (temperatures) at time t m .The main idea of the proper orthogonal decomposition is to obtain the most typical or characteristic structure /(x) among these T m (x) "m.The maximization of leads to which can be rewritten in the form Defining the vector / such that its i-component is /(x i ), Eq. ( 7) results in the eigenvalue problem (8), whose eigenvectors related to the highest eigenvalues define the characteristic structure of T m (x) where the two points correlation matrix is given by which is symmetric and positive definite.If we define the matrix Q containing the discrete field history then it is easy to verify that the matrix c in Eq. (8) results

A posteriori reduced modelling
We solve the eigenvalue problem defined by Eq. (8) selecting the n eigenfunctions / k associated with the eigenvalues belonging to the interval defined by the highest eigenvalue and that value divided by a large enough value (10 8 in our simulations).In practice n is much lower than N. Now, we could try to use these n eigenfunctions / k for approximating the solution T m (x) "m.For this purpose we need to define the matrix Now, if we consider the linear system of equations resulting from the semi-implicit discretization of the thermal model defined by Eq. ( 4) expressing Eq. ( 13) results and by multiplying both terms by B T we obtain where the size of When n ( N, as is the case in numerous physical models, the solution of Eq. ( 16) is preferred because its reduced size.
Remark 1. Eq. ( 16) can be also deduced by introducing the approximation ( 14) into the Galerkin variational form related to the partial differential equation.If trial and test functions in Eq. (2) are approximated using the form (14) then Eq. ( 16) is easily deduced.
Remark 2. An alternative technique to reduce the size of the eigenvalue problem lies in the application of the snapshot proper orthogonal decomposition.This technique is based on finding the significant modes from the application of the POD, but assuming that those modes can be written as a linear combination of the M snapshots that were used to define the decomposition.The main advantage of this strategy is that theses modes result from the eigenproblem defined by Only the eigenvectors related to large enough eigenvalues are retained.From these eigenvectors the reduced approximation functions are computed using the fact that these functions are linear combination of the snapshots, i.e.
Remark 3. In the previous analysis the reduced basis was built from the computed unknown field evolution that was carried out solving the discrete evolution problem (4).Thus, one could ask about the interest of a such approach.There are two kind of approaches widely considered.The first approach consists in solving the non-reduced model (Eq.( 4)) in a short time interval, allowing the extraction of the characteristic functions and then the definition of the reduced approximation basis that is then used to perform the solution of the reduced evolution model (Eq.( 16)) in larger time intervals with the associated computing time savings.The other approach consists in solving the nonreduced model (Eq.( 4)) in the whole time interval, whose solution allows defining the reduced approximation basis that then could be used for solving ''similar'' models, as the ones that involve for example slight variations in some material parameters or in the boundary conditions.Some recent advances in such approaches can be found in [16,15,18,5,7,12] and the references therein.

Enriching the approximation basis
Obviously, accurate simulations require an error evaluation as well as the possibility of adapting the approximation basis by introducing new functions able to describe the solution features.Ryckelynck proposed in [17] an adaptive procedure, able to construct or enrich the reduced approximation basis.For this purpose, he proposed to enrich the reduced approximation basis by adding some Krylov subspaces generated by the equation residual.Despite the fact that this enrichment tends to increases the number of approximation functions, when it is combined with a POD decomposition that continuously reduces this number, the size of the problems is quickly stabilized.This strategy, which was successfully used in [1] as well as in [19], is summarized in the present section.
We assume the solution accurately described in the interval ]0, t s ] by using the reduced basis B. Now, the solution evolution is performed in the time interval ]t s , t s+1 ] solving Eq. ( 16) making use of the reduced approximation basis defined by matrix B: When time t s+1 is reached, a control step is performed in order to evaluate the accuracy of the solution computed using the reduced basis.The control step is performed only at the end of each time interval.We assume that t s+1 À t s = M • Dt and consequently the residual at t s+1 , R M , can be computed from If the norm of the residual is small enough kR M k < kBf M k (being a small enough parameter) the computed solution can be assumed as good, and the time integration goes on in ]t s+1 , t s+2 ] using Eq. ( 19) without changing the reduced approximation basis.On the contrary, if kR M k P kBf M k, the approximation must be improved.For this purpose, we propose to enrich the reduced approximation basis by introducing the just computed residual (or some Krylov's subspaces generated by the residual, as proposed in [17]): and now, the evolution is recomputed in the ]t s , t s+1 ] using Eq. ( 19) with the just updated reduced basis B. When the convergence criterion is satisfied (i.e.kR M k < kB f M k) we apply a POD on the entire past history ]0, t s+1 ] for defining the reduced basis able to represent all the past evolution of the unknown field, which results in an updated B that can be considered as the optimal reduced basis to describe the evolution of the unknown field in ]0, t s+1 ].Using the just updated reduced basis B the time evolution given by Eq. ( 19) is performed in ]t s+1 , t s+2 ] and the solution accuracy is checked at t = t s+2 .If a basis enrichment was needed at the end of the previous time interval, at t = t s+1 , then the length of the present time interval is reduced according to t s+2 À t s+1 = (t s+1 À t s )/2, and if not enrichment was needed, then the time interval length is increased according to t s+2 À t s+1 = 2(t s+1 À t s ) (see [1] for more details on the time interval length adaptation).
The enrichment tends to increase the size of the reduced basis, but the POD reduces its size.The combination of both procedures allows to stabilize the size of the model as we noticed in some of our former works [1].

Parallel time integration
The reduction technique just described can be coupled with parallel time integration, naturally addressed by the too small time steps required in the simulation.For performing parallel time integration the entire time domain is partitioned in a number of time intervals.
A well established technique for parallel time integration, the ''parare ´elle'' approach, was proposed in [14].In this approach the solution in each time interval is computed from a trial initial condition.The solution jumps at the time intervals boundaries are used to update the prescribed initial conditions related to each time interval.The iteration procedure stops when no solution discontinuities are noticed across the time intervals boundaries.As we can notice this strategy needs for an iteration procedure and adequate criteria to update the initial conditions at each time interval as well as to evaluate the convergence.
In that follows we propose an alternative approach, slightly different, based on the use of model reduction and some results of ODE theory.When we consider a reduced modelling the computing time related to both approaches (our approach and the ''parare ´elle'' one) is similar.However, our approach does not involve neither an iteration scheme nor the necessity of introducing convergence criteria.

Solving systems of linear ordinary differential equations
In this section we consider the general discrete form of a linear PDE from which we can identify the solutions of the homogeneous T h and complete T c problems: The general solution of Eq. ( 22) could be computed from the general solution of the homogeneous problem and a particular solution of the complete one, that is The general solution of the homogeneous problem T h (t) results where T i h ðtÞ is computed by solving the homogeneous system related to the initial condition where the jth component of vector d i results the Kroenecker's delta d ij .
On the other hand, the particular solution T c (t) can be computed by integrating the second differential equation in (23) from any initial condition, for example T c (t = 0) = 0. Now, the general solution is perfectly defined where the coefficients a i must be computed in order to satisfy the initial condition: At present, this approach, originally proposed in [10], does not seem very useful, but we will prove in the following section that this approach allows considering parallel time integration.

Introducing reduced representations
The just described integration scheme can be generalized for solving the evolution problem defined in the reduced approximation basis.For this purpose we consider instead Eq. ( 22) its reduced counterpart: From a similar reasoning, we can write where functions f i h ðtÞ are solution of with f c (t = 0) = 0.
The n coefficients b i (related to the n functions defining the reduced approximation basis) must be calculated to verify the initial condition Now, two possibilities exist.The first one is based on the weak imposition of the initial condition whereas the second one is based on adding the initial condition to the reduced approximation basis, that is which allows to write from Eq. ( 34): b T = (1, 0, . . ., 0); from which it can be noticed that only the solution related to f 1 h ðtÞ is required.

Time interval partitioning in the linear case
Now, we consider the time interval [0, t max ] partitioned into P intervals where t 0 = 0 and t P = t max .We firstly assume that all the time intervals have the same reduced approximation basis B. The time integration can be performed simultaneously in the different time intervals according to the procedure described in the previous section.
If the reduced approximation basis contains n functions, then n + 1 problems must be integrated in each time interval, n related to the general solution of the homogeneous system and the last one associated with a particular solution of the complete system.Thus, for each time interval We must notice that the resolution of Eqs. ( 38) and (39) can be performed simultaneously (in a parallel computing platform) for the different time intervals I p .Now, for computing all the coefficients b p i we impose the continuity conditions between the different time intervals as well as the initial condition that implies as discussed at the end of the previous section either or if the first function of the reduced approximation is the initial condition.It must be noticed that when the general solution of the homogeneous equation and a particular solution of the complete one have been defined on each time interval, then Eq. (44) (or Eq. ( 45)) combined with Eq. (43) allow to define immediately the numerical solution on the entire time domain.We must also recall that the computation of the solution in each time interval only involves the resolution of n linear transient models of size n (n being the number of degrees of freedom, that is, the number of functions defining the reduced approximation basis, which usually is of the order of tens) and that solution can be performed in parallel on each time interval.
If we are using the same reduced approximation basis in all the time intervals (considered with the same length), and the same number of time intervals are addressed to the different computer processors, the charge of each processor will be the same.

Enriching the approximation bases
As previously discussed the approximation basis can be improved by adding the residual or some Krylov subspaces related to the governing equations residual.Thus, with the solution determined according to Eq. (40) in each time interval I p , the residual R p (t = t p ), can be computed "p P 1 from

Now two simple possibilities exist
(1) To perform a global enrichment and recomputed all the solutions until convergence, that is, while (2) To perform a local enrichment "p from which the different solutions must be recomputed in each time interval, but now matrix and vectors in Eqs. ( 38) and (39) depend on the time interval considered from the dependence of these entities on the reduced approximation bases according to Eq. (30).
In the case of considering different reduced approximation basis in each time interval, the transmission conditions between neighbor intervals must be reformulated.There are different possibilities of performing this transmission, the simplest one based on enforcing "p 2 [1, P À 1]: that must be satisfied in a weak sense In order to keep the reduced approximation order, after some bases enrichment we can perform a POD decomposition in each interval to extract the significant information on the solution evolution.The eigenfunctions selection criterion is based in the relative value of the associated eigenvalues (with respect to the highest eigenvalue) as described in Section 1.

Time interval partitioning in the non-linear case
In the non-linear case (out of the scope of the present work due to the negligible influence of non-linearities in the thermal processes here considered) the procedure described in Section 2.1 cannot be applied.However, as proposed in [9] it could be applied if the solution of that non-linear problem is searched from the solution of a sequence of linear problems using standard linearization strategies (fixed point, Newton).In what follows the simplest scheme which is based on the fixed point linearization, is described and used.
Thus, coming back to Section 2.3, and assuming the same reduced approximation basis for all the intervals I p , we must compute, for each time interval I p , at iteration l P 1 In the previous expressions we have assumed a dependency of the thermal conductivity and the specific heat on temperature.The solution in each interval, after the imposition of the transmission conditions, could be defined either by using that can be integrated using an explicit scheme from the initial condition expressed as We must notice that, as in the linear case, the resolution of Eqs. ( 52) and (53) can be performed simultaneously (in a parallel computing platform) for the different time intervals I p .The iteration procedure stops when the fixed point is reached, that is kT p;ðlþ1Þ ðtÞ À T p;ðlÞ ðtÞk dt < : If one consider the solution of Eqs. ( 52) and (53) with the solution reconstruction given by Eq. ( 55) the proposed strategy always converge, at most in P iterations, because in the first iteration the solution in the first interval is perfectly determined from the initial condition, which implies exact initial condition for the second interval at the second iteration, and so on.In general the fixed point algorithm converges quickly, and in this case if the number of time intervals P is large enough, significant CPU times savings could be attained.However, if the number of time interval is lower than the number of fixed point iterations needed for reaching convergence, then the efficiency of the method becomes similar to the efficiency of a sequential solution.
The number of iterations needed to reach the convergence depends on the treated case, however, our present numerical experiences as well as the ones reported in [9] seem indicating that this number is around 5 (for moderate non-linearities).

Numerical example
For illustrating both procedures just described (the one based on the model reduction and the one based on the time partitioning) we consider in this section a simple example concerning a 2D square domain X=]0, 1[ • ]0, 1[ and the time interval [0, t max = 30].We consider x = (x, y), o T X = {x:y = 1} and o Q X the complementary part of the domain boundary.The temperature is then suddenly prescribed on o T X "t > 0, according to that induces a thermal chock at t = 0.
In o Q X we assume a null heat flux across the boundary.The initial condition imposed at the initial time is given by T(x, t = 0) = T 0 (x) = 1.We also assume that the medium is homogeneous and isotropic, which reduces the thermal conductivity tensor to a scalar.In the linear case we assume constant conductivity K = Id (Id being the unit tensor) and specific heat qC = 1.
The domain mesh consists of 20 • 20 four nodes quadrilateral elements, the time step is fixed to Dt = 0.1 and a fully implicit strategy is considered in the time integration.In Fig. 1 we depict the temperature evolution at four different times.During the time evolution, the history matrix Q (Eq.( 10)) is built-up taking the temperature each ten time steps.Thus the resulting Q matrix contains 27 temperature snapshots.Now, the POD decomposition can be applied, solving the eigenvalue problem defined by Eq. ( 8), from which it results only 4 eigenfunctions related to the 4 most significant eigenvalues (the remaining eigenvalues being lower than the highest one multiplied by 10 À8 we decided to neglect the influence of the associated eigenfunctions in the solution evolution approximation).
For the sake of simplicity in the prescription of the initial condition we propose to add to this four functions the one related to the initial temperature according to Eq. (36).In Fig. 2 we depict the 3 first normalized eigenfunctions, as well as the one corresponding to the initial condition.If the solution evolution is now recomputed using the reduced approximation basis that only involves 5 degrees of freedom, the computed solution is in perfect agreement with the one computed using the finite element approximation basis.For proving it, we depict in Fig. 3 the comparison between T FEM (x, t = t max ), and the one computed in the reduced approximation basis T RED (x, t = t max ).As it can be noticed the maximum differences are of order 10 À9 .
If the time domain is partitioned in several intervals the computed solution is, as expected, exactly the same than the one computed using a single time interval, because both computations have been carried out using the same reduced basis.Now, in order to validate the adaptation scheme we consider the reduced basis just computed applied to the solution of a different problem.The only differences between the problem that served to compute B and the present one concerns the essential boundary condition that at present is given by We are solving the problem using two time intervals with the same length.The computed solution at t max using the reduced approximation basis previously computed involves a slight deviation with respect to the one computed by using the FEM.Fig. 4 compares both solutions: T FEM (x, t max )  and T RED (x, t max ).The error norm considering the natural norm L 2 (0, t max ; L 2 (X)) results kT RED À T FEM k L 2 % 0:1.The reduced basis used for computing the temperature evolution was associated with the boundary condition given by Eq. ( 58) that does not correspond with the one prescribed at present (Eq.( 59)), which justifies the slight deviation noticed.Now, in order to reduce this deviation we proceed to enrich the reduced approximation basis from the residuals computed at the end of each time interval according to Eq. ( 47) Now, the entire whole evolution is recomputed and the basis enriched until convergence.
The convergence is attained in 3 iterations, which implies the introduction of 6 new approximations functions.The comparison of the FEM and the reduced basis solutions at t max is depicted in Fig. 5.The associated error norm is reduced in the order of 100.

High resolution homogenization
Due to the microscopic heterogeneity, the macroscopic thermal modelling needs an homogenized thermal conductivity which depends on the microscopic details (fibers arrangement in the microscopic representative volume related to each point in the macroscopic domain).
To compute this homogenized thermal conductivity in the linear case, one could isolate a representative volume element X rve (we do not consider here the problematic related to the definition of such representative volume which is addressed for example in [13]) and assume that the microstructure is perfectly defined at such scale.Thus, the microscopic conductivity is assumed known at each point in the microscopic domain k(x).
We can define de macroscopic temperature gradient g Grad T ¼ hGrad T i as We also assume the existence of a localization tensor L(x) such that Grad T ðxÞ ¼ LðxÞ g Grad T : Now, we consider the microscopic heat flux q according to the Fourier's law from which the macroscopic counterpart Q results Fig. 3. Finite element and reduced approximation basis solution at t = t max : (a) Computed solution using the reduced approximation basis; (b) T FEM (x, t = t max ) À T RED (x, t = t max ).from which the homogenized thermal conductivity can be defined as As k(x) is perfectly known everywhere in the representative volume element, the definition of the homogenized thermal conductivity tensor just requires the computation of the localization tensor L(x).For this purpose we consider in the 3D case the solution of the three boundary value problems related to the steady state heat transfer model in the microscopic domain, defined by DivðkðxÞGrad T 1 ðxÞÞ ¼ 0; DivðkðxÞGrad T 3 ðxÞÞ ¼ 0 It is easy to prove [11] that these three solutions verify where (AE) T denotes the transpose.Thus, the localization tensor results finally From the computational point of view the main difficulty concerns the solution of the three boundary value problems using very fine meshes required to represent accurately all the microscopic details.

Separated representation and tensor product approximation basis
We first consider a representative volume element that for composite materials applications must contain numer-ous fibers.Due to the characteristic size of those fibers a very fine mesh will be required for solving accurately the thermal models related to Eqs. ( 66)-( 68).Thus, it is usual to define 3D models consisting of thousands of millions of nodes.Nowadays, the performance of computers allow to solve, in a reasonable time, problems rarely exceeding several millions of degrees of freedom.
In order to circumvent this difficulty we propose to use a separated representation in combination with a tensor product reduced approximation basis, which allows to a discrete model of size N • D (N being the number of nodes in each space direction and D the dimension of the space) instead the N D characteristic of grid-based discrete models.Thus, if we suppose N % 10 5 and D = 3, the size of the discrete model will be 3 • 10 5 instead 10 15 .
In the following paragraphs we summarize the main elements of this technique that we have successfully applied for solving some multidimensional models usually encountered in the kinetic theory models of complex fluids in steady [2] and transient [3] regimes.
We consider the steady conduction heat transfer problem in X=]0, L[ 3 , the temperature vanishing on the domain boundary o T X oX, T(x 2 o T X) = 0 (the non-homogeneous case will be addressed later) where the temperature T and the microscopic thermal conductivity tensor k depend on the space coordinates x = (x 1 , x 2 , x 3 ).The problem solution can be assumed in the form where F kj 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 of approximation functions, i.e.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 grid).At the present stage of development we consider that the 1D meshes are fine enough for capturing all the solution details.The use of adaptive separated representations in the context of multiresolution analysis (e.g. using wavelets) able to capture solution details and to account for different characteristic length scales, as used in the sparse grid framework [6], constitutes a work in progress.In any case, as proved in [2,3], this kind of approximation can describe accurately solutions that exhibit boundary layers.Now, an appropriate numerical procedure is needed for computing the coefficients a j as well as the J 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 kj ("j 2 [1, . . ., n]; "k 2 [1, . . ., 3]) known (verifying the boundary conditions), the coefficients a j can be computed by introducing the approximation of T into the Galerkin variational formulation associated with Eq. ( 71): Thus, using the approximation of T T ðx 1 ; x 2 ; and Now, we assume that f(x) and k(x) can be written in the form This kind of separated representation can be easily performed using singular value decomposition (in 2D) or using an alternating least squares technique in the general multidimensional case [4].For this purpose we only need to know the values of the scalar or tensorial field in a cloud of nodes.Eq. ( 77) involves integrals of a product of 3 functions each one defined in a different dimension.Let Q 3 k¼1 g k ðx k Þ be one of these functions to be integrated.The integral over X can be performed by integrating each function in its 1D definition interval and then multiplying the 3 computed integrals according to which makes possible a fast numerical integration.Now, due to the arbitrariness of the coefficients a Ã j , Eq. (77) allows to compute the n coefficients a j , solving the resulting linear system of size n • n.This problem is linear and moreover n rarely exceeds the order of tens.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 T at iteration n given by Eq. ( 75) we compute the residual R related to Eq. ( 71) If R < ( being a small enough parameter) the iteration process stops, yielding the solution T(x) given by Eq. (75).Otherwise, the iteration procedure continues.The integral in Eq. ( 80) 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 a j just computed the approximation basis can be enriched by adding the new functions Q 3 k¼1 F kðnþ1Þ ðx k Þ.For this purpose we solve the non-linear Galerkin variational formulation related to Eq. (71): using the approximation of T given by The solution of Eq. ( 81) can be expressed from the stationarity of functional J(T) that corresponds to Eq. ( 81) by putting dT = T * .Now, taking the variation of T according to its expression given by Eq. ( 82), where the variation of the known functions F kj vanishes, we have that can be rewritten as This leads to a non-linear variational problem, whose solution allows to compute the 3 functions R k (x k ).Functions F k(n+1) (x k ) are finally obtained by normalizing, after convergence of the non-linear problem, the functions R 1 , R 2 , R 3 .
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 interpolation.If we assume that p k nodes are used to define the interpolation of function R k (x k ) in the interval [0, L], then the size of the resulting discrete non-linear problem is P k¼3 k¼1 p k .The price to pay for avoiding a whole mesh in the domain is the solution of some non-linear problems of size P k¼3 k¼1 p k .However, the size of those non-linear problems remains moderate and no particular difficulties have been found in its solution.
Different non-linear solvers were analyzed in [2,3].The first strategy was based on the Newton scheme and the second in an alternating directions strategy.The numerical results presented in this paper have been computed using the alternating directions method.It is important to note that the discrete problems resulting at this enrichment stage when one is using an alternating directions procedure are three-diagonal.Remark 4. Steps 2 and 3 constitute a convergent procedure but we have noticed experimentally that the consideration of step 1 accelerate the convergence in numerous cases, as those concerning kinetic theory models widely considered in our former works [2,3].By this reason we decided to keep this projection step in the global algorithm, without a significant incidence on the computing time.

Numerical example
We consider the steady conduction heat transfer problem defined in a 2D square domain X = ]0, L[ 2 where the conductivity is assumed isotropic and where the tempera-ture is prescribed on its boundary o T X oX according to T(x 2 o T X) = T g (x) = x 1 (which corresponds to problem (66)) where the temperature T and the microscopic thermal conductivity k depend on the space coordinates x = (x 1 , x 2 ).The associated variational formulation results The problem solution can be written in the form Now, assuming that the conductivity can be separated (using singular value decomposition or an alternating moving least squares technique): the integrals related to the 3 steps of the procedure described in the previous section can be performed, and after convergence the solution T(x) is obtained from e T ðxÞ according to Eq. (89).
Fig. 6 illustrate the conductivity at each point x 2 X as well as the different functions k 1h (x 1 ) and k 2h (x 2 ) involved in its separated representation Finally, Fig. 7 depicts the computed temperature fields e T ðxÞ and T(x), respectively, where x 1 x, x 2 y, F 1j F j and F 2j G j .In this calculation 10 4 nodes were considered in each direction.Thus, the size of the problem was 2 • 10 4 instead the 10 8 that should been required in the context of the finite element method for computing a solution with an equivalent accuracy (we noticed in [3] that this strategy exhibits a fourth order convergence rate).The separated representation only involves 5 one-dimensional functions, all them represented in Fig. 7.

Conclusion
This paper explores the ability of model reduction coupled with parallel time integration for solving thermal models.The main interest of that approach lies in the fact that when the domain mesh involves tremendous number of degrees of freedom being the time step also reduced because the existence of fast evolving boundary conditions, the combination of both strategies (the first one that allows reduce the number of degrees of freedom and the second one making possible the use of parallel computer platforms for the time integration) seems to be an appealing choice.
Moreover, the use of separated representation and the definition of tensor product reduced approximation basis allows to employ high resolution homogenization techniques, because in this case, even very fine 3D models result perfectly tractable.

Fig. 4 .
Fig. 4. Finite element and reduced approximation basis solution at t = t max before the basis adaptation: (a) Computed solution using the reduced approximation basis; (b) T FEM (x, t = t max ) À T RED (x, t = t max ).

Fig. 5 .
Fig. 5. Finite element and reduced approximation basis solution at t = t max after the basis enrichment: (a) Computed solution using the reduced approximation basis; (b) T FEM (x, t = t max ) À T RED (x, t = t max ).