J.-C Passieux
email: passieux@lmt.ens-cachan.fr

P Ladevèze
email: ladeveze@lmt.ens-cachan.fr

D Néron
email: neron@lmt.ens-cachan.fr

A scalable time-space multiscale domain decomposition method: adaptive time scale separation

Keywords: scalability, multiscale in time and space, domain decomposition, parallel computing, model reduction

This paper deals with the scalability of a time-space multiscale domain decomposition method in the framework of time-dependent nonlinear problems. The strategy which is being studied is the multiscale LATIN method, whose scalability was shown in previous works when the distinction between macro and micro parts is made on the spatial level alone. The objective of this work is to propose an explanation of the loss-of-scalability phenomenon, along with a remedy which guarantees full scalability provided a suitable macro time part is chosen. This technique, which is quite general, is based on an adaptive separation of scales which is achieved by adding the most relevant functions to the temporal macrobasis automatically. When this method is used, the numerical scalability of the strategy is confirmed by the examples presented.

Introduction

In computational mechanics, the simulation of the behavior of complex structures in which two or more very different scales can be present in both space and time is a challenging question. A typical engineering example is that of a relatively large structure in which local cracking or local buckling occurs [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods in structural mechanics[END_REF]. Another typical engineering problem is related to today's deep interest in material models described on a scale smaller than that of the macroscopic structural level, such as composite materials [START_REF] Ladevèze | Computational damage micromodel of laminated composites[END_REF]. In such situations, the local solution involves short-wavelength phenomena in both space and time which require the use of complex models describing the material on a very refined scale. As a result, classical finite element codes lead to systems with very large numbers of degrees of freedom whose resolution would generally be excessively expensive. Therefore, one of today's main challenges is to derive computational strategies capable of solving such engineering problems through true interaction among the scales, both in space and in time.

As far as space is concerned, one of the earliest strategies consisted in applying the theory of the homogenization of periodic media initiated by Sanchez-Palencia [START_REF] Sanchez-Palencia | Non homogeneous media and vibration theory[END_REF]. Similar computational approaches can be found in [START_REF] Devries | Homogenization and damage for composite structures[END_REF][START_REF] Fish | Computational plasticity for composite structures based on mathematical homogenization: Theory and practice[END_REF][START_REF] Oden | Hierarchical modeling of heterogeneous solids[END_REF][START_REF] Feyel | A multilevel finite element (FE2) to describe the response of highly non-linear structures using generalized continua[END_REF][START_REF] Temizer | An adaptive method for homogenization in orthotropic nonlinear elasticity[END_REF]. First, the macroproblem leads to effective values of the unknowns; then, the microsolution must be calculated locally using localization operators. The fundamental assumption, besides periodicity-which is not required, thanks to thermodynamical arguments such as Hill-Mandel's criterion-, is that the ratio of the two scales is small. The boundary zones, in which the material cannot be homogenized, can be treated using specific techniques [START_REF] Devries | Homogenization and damage for composite structures[END_REF]. Another type of approach uses the macroscale as the reference scale, and requires the microscale to enrich the model only where the macroscale is not sufficient. The microscale may be associated with a more refined description [START_REF] Ibrahimbegović | Strong coupling methods in multi-phase and multi-scale modeling of inelastic behavior of heterogeneous structures[END_REF], a different model [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF], or even an analytical approach [START_REF] Hughes | The variarional multiscale-a paradigm for computational mechanics[END_REF]. Space can also be enriched using some a priori knowledge of the solution [START_REF] Melenk | The partition of unity finite element method: basic theory and applications[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Strouboulis | The generalized finite element method[END_REF]. The last family of methods, which includes domain decomposition methods [START_REF] Mandel | Balancing domain decomposition[END_REF][START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF][START_REF] Ladevèze | A micro/macro approch for parallel computing of heterogeneous structures[END_REF][START_REF] Alart | A scalable multiscale LATIN method adapted to nonsmooth discrete media[END_REF][START_REF] Pebrel | A nonlinear dual domain decomposition method: application to structural problems with damage[END_REF] and multigrid methods [START_REF] Briggs | A multigrid tutorial[END_REF][START_REF] Fish | Multigrid method for periodic heterogeneous media II. multiscale modeling and quality control in multidimensional case[END_REF], starts from the microscale and uses the larger scales to accelerate the convergence of an iterative algorithm.

Among the relatively few works which have been devoted to multi-time-scale computational strategies, one can find multi-time-step methods [START_REF] Gravouil | Multi-time-step and two-scale domain decomposition method for nonlinear structural dynamics[END_REF][START_REF] Combescure | A numerical scheme to couple subdomains with different timesteps for predominantly linear transient analysis[END_REF] and parallel time integrators [START_REF] Belytschko | Stability of multi-time step partitioned integrators for firstorder finite element systems[END_REF][START_REF] Bottasso | Multiscale temporal integration[END_REF][START_REF] Farhat | Time-decomposed parallel time-integrators: theory and feasibility studies for fluid, structure, and fluid-structure applications[END_REF], which deal with different time discretizations and integration schemes. For multiphysics problems, coupling among time grids has been proposed in [START_REF] Néron | A computational strategy for poroelastic problems with a time interface between coupling physics[END_REF]. In these approaches, each region of space or each physics is described using a single scale. Other approaches involve a true two-scale description of time phenomena through the introduction of "micro-macro projectors" [START_REF] Maday | A parareal in time procedure for the control of partial differential equations[END_REF]. Time homogenization methods have been introduced for dealing with periodic loading histories [START_REF] Guennouni | On a computational method for cycling loading: the time homogenization[END_REF][START_REF] Cognard | A large time increment approach for cyclic plasticity[END_REF] and time-space homogenization for analyzing multiple physical processes interacting at multiple temporal and spatial scales [START_REF] Yu | Multiscale asymptotic homogenization for multiphysics problems with multiple spatial and temporal scales: a coupled thermoviscoelastic example problem[END_REF].

This paper focuses on the separability of scales, especially time scales. The method which serves as the starting point of this discussion is a mixed, multilevel domain decomposition method based on the LATIN method [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF]. This method leads to resolutions both on the refined scale (the microscale) and on the coarse scale (the macroscale) with respect to time as well as space. The problem on the microscale is solved using the iterative LATIN solver while the purpose of the macroresolutions is to accelerate the convergence of this iterative algorithm. The choice of the definition of the macro description has no influence on the solution after convergence, but affects the convergence rate alone. Besides, different time discretizations can be used from a subdomain to another.

Regarding the space variable, the choice of the macrospace proposed in [START_REF] Ladevèze | A micro/macro approch for parallel computing of heterogeneous structures[END_REF] was found to be optimal from a numerical scalability point of view. Indeed, the spatial macrobasis contains the resultants and moments of the connecting forces at the interface, and due to Saint Venant's principle the microcomplement (whose resultant and moment are zero) has only a local incidence in the few substructures surrounding that interface. Thus, the spatial macrospace, in spite of being associated with a very small number of unknowns at each interface, has a strong physical meaning and leads to a numerically scalable algorithm (as long as, regarding time, the scales reduce to the microscale).

Regarding time, the natural choice which consists in defining a priori the macroscale as the slow scale and the microscale as the rapidly evolving scale in time fails to preserve that optimality, which makes the method nonscalable in the general case. A seemingly natural improvement would be to enrich the macrobasis using an h or p method, but we will see on some examples that that is not sufficient to restore the numerical scalability property.

The objective of the paper is not to compare the LATIN method to other multiscale strategies, but to find a remedy for this scalability problem. In order to do that, we propose to introduce an adaptive approach to time scale separation which consists in enriching the temporal macrobasis using the most suitable modes. The idea is to view the macro time part of a quantity as its projection onto a reduced basis whose vectors are not necessarily known at the beginning of the calculation. Thus, the basis is generated automatically. It depends on the problem and on the loading, and it can also evolve during the iterative resolution of a given problem.

Indeed, an automatic enrichment technique enables one to update this basis during the calculation. This enrichment technique consists, at a particular iteration of the algorithm, in selecting a set of microresiduals whose resultants and moments are nonzero. According to Saint Venant's principle, these functions are representative of microphenomena which have global influence, but are not taken into account by the macroproblem. All these functions could a priori be added to the temporal macrobasis, but this would lead to prohibitive computation costs. A technique based on the Proper Orthogonal Decomposition (POD) method ( [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF]) is used to extract from these functions the most representative elements in order to add them to the temporal macrobasis. Thus, the space of the macro time fields is improved and adapted to the problem as the iterations go on.

The article is structured as follows. In Section 2, the reference problem is introduced along with a brief review of the main aspects of the multiscale method which serves as the reference for the work presented here. In Section 3, the choice of the spatial macrobasis is reviewed and the numerical scalability issue is explained and studied through an example. In Section 4, we show that in some cases the natural choice of the macro time space is not optimum and that h-adaptivity or padaptivity of the macrospace fails to solve that problem, which shows that, in the general case the method is not numerically scalable. In Section 5, we introduce our proposed adaptive technique in order to update the temporal macrobasis. We also present the selection technique which ensures that the additional cost associated with our method remains small. Finally, in Section 6, the numerical scalability of the adaptive method is illustrated using a 3D viscoelastic problem with frictional contact.

The time-space multiscale LATIN method

This section gives a brief review of the main aspects of the multiscale computational strategy. Further details can be found in [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF][START_REF] Ladevèze | On multiscale computational mechanics with time-space homogenization[END_REF].

The reference problem

Let us consider the evolution of a structure defined over a time-space domain I ×Ω, where I = [0, T ]. The structure is subjected to prescribed body forces f d , traction forces F d , and prescribed displacements U d .

The structure is viewed as an assembly of substructures and interfaces. Similarly, the time domain I is divided into a small number of coarse intervals

I C i = [t C i , t C i+1 ]
, each in turn being divided into more refined subintervals

I f j = [t f j , t f j+1 ].
Let Φ EE ′ denote the interface between substructures Ω E and Ω E ′1 . This interface is characterized by the restriction to Φ EE ′ of the displacement fields (W E , W E ′ ) and the force fields

(F E , F E ′ ).
The behavior of the interface is characterized by the introduction of a relation among these quantities which is detailed in [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF]. (For frictional contact, see [START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF].)

Two-scale description of the unknowns

The scale separation between a micro part m and a macro part M takes place only at the interfaces and, ∀F M ⋆ ∈ F M E , is defined by:

I C i ×∂Ω E ( Ẇ M E -Ẇ E ) • F M ⋆ dSdt = 0 and Ẇ m E = Ẇ E - Ẇ M E (1)
The spaces F M E and W M E can be chosen arbitrarily. This choice, which will be discussed in Sections 3 and 4, has a strong influence on the scalability of the method.

A major point of the strategy, which grants it its multiscale character, is that the set of the macro forces F M = (F M E ) Ω E ⊂Ω is required to verify the transmission conditions at the interfaces a priori at each iteration.

The algorithm

This problem is solved using the LATIN method [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF], a general iterative nonlinear resolution technique for time-dependent problems which operates globally over the entire time-space domain. One iteration consists of two stages, called the "local stage" and the "linear stage", in which solutions verifying the nonlinear constitutive relation (defined by a space Γ) and a group of admissibility equations called A d are built alternatively. In order to close the problem, one needs to introduce what we call the "search directions" E + and E -, which are detailed in the above references.

The local stage consists in a set of nonlinear problems which are solved independently at each point of the discretization of each subdomain Ω E and each interface Φ EE ′ , and for each integration point of the time domain I. The resolution of that stage is straightforward and leads to a solution ŝ = {ŝ E } E .

In the linear stage, one seeks a solution which verifies the admissibility conditions over each Ω E by following the search direction, which involves an operator which can be interpreted as a linearized constitutive relation. Therefore, the problem is linear, but not independent within each substructure because the admissibility of the macroforces couples all the subdomains. This problem leads to a solution s = {s E } E . The admissibility conditions are written in the weak sense, using

Lagrange multiplier ˙ W M E whose admissibility is expressed, ∀ ˙ W M ⋆ ∈ W M ⋆ E , by
:

Ω E ⊂Ω I C i ×∂Ω E ˙ W M ⋆ E • F E dSdt I C i ×∂Ω E ∩∂ 2 Ω - ˙ W M ⋆ E • F d dSdt = 0 (2)
For the sake of simplicity, let us introduce the linear operator L E , associated with the resolution of the microproblem, which maps its right-hand side to solution s E , e.g.:

s E = L E (ŝ E + ˙ W M E ) (3)
At this stage, ˙ W M E is unknown, but since this microproblem is linear its resolution can be divided into two parts. First, one can calculate s 1 E , the solution of the microproblem associated with the known solution of the previous local stage ŝE . The microproblems in this first set can be calculated independently.

s 1 E = L E (ŝ E ) (4)
Knowing s 1 E , one can calculate the macro part of the corresponding force distribution F M,1 . We will see that this quantity is used to assemble the right-hand side of the macroproblem. The remainder of the microproblem is:

s 2 E = L E ( ˙ W M E ) (5)
If solution s 2 E is known (which is not yet the case), one can calculate the associated macro force distribution F M,2 . In particular, Problem (5) can be written for the macro part of the forces alone:

F M,2 = L F E ( ˙ W M E ) (6)
Since Equation ( 6) maps a small-dimension space to the same space of macro fields, operator L F E is discrete by nature, L F E , and can be calculated explicitly at reasonable cost.

F M,2 = L F E ˙ W M E (7)
Operator L F E , called the homogenized operator over the time-space subdomain Ω E × I C i , can be viewed as a time-space macro/mixed Schur complement of subdomain Ω E × I C i . Its calculation requires a series of resolutions of Problem [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF] in which ˙ W M E takes the values of each basis vector of W M E successively. This operator depends only on the choice of the macrobases and on the parameters of the search directions. Therefore, it can be constant over a large number of iterations provided that the search directions remain unchanged.

Now the solution of the linear stage as a function of the macroforces is:

F M = F M,1 + F M,2 = F M,1 + L F E ˙ W M E (8)
Finally, the macroproblem is obtained by introducing Form (8) of the macro forces field into Relation (2). This linear time-space problem is defined over the whole set of interfaces and the entire coarse subinterval I C i . One can prove that this problem has a unique solution.

Ω E × I C i , ∀E. • Convergence criterion η = s -ŝ end Algorithm 1:
The multiscale LATIN method Algorithm 1 summarizes the organization of the different stages of the method. For further details, see [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF][START_REF] Ladevèze | On multiscale computational mechanics with time-space homogenization[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF].

The efficiency of this method in terms of computing time compared with classical approaches is discussed in [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF][START_REF] Ladevèze | A micromacro and parallel computational strategy for highly heterogeneous structures[END_REF].

Another important aspect of the LATIN method which is not presented here is the use of Proper Generalized Decomposition (PGD) methods ( [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF][START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF]) which enable a drastic reduction of the amount of data to be handled and the number of calculations. For further information, the reader can refer to [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF].

Choice of the spatial macrobasis

Scalability

Domain decomposition methods were initially developed in order to take advantage of parallel computer architectures. The idea is simple: if a single processor is capable of solving a problem of up to size N , using n processors and an appropri-ate method, one can hope to solve a problem of size N × n.

The speedup, denoted s(n), is defined as the ratio of the sequential execution time to the parallel execution time. The efficiency, denoted e(n), is the ratio of the speedup to the number of processors n: e(n) = s(n) n An algorithm is considered to be efficient if e(n) is close to 1. If the efficiency of an algorithm is stable when n increases, that algorithm is said to be numerically scalable. Scalability characterizes the robustness of the approach when the number of processors (i.e. the number of substructures) increases. In order to deal with large problems, it is absolutely necessary to use a scalable method. In the context of domain decomposition methods, scalability is linked to the dependence of the convergence rate on the number of substructures. This paper deals with the influence of the choice of the spatial and temporal macrospaces (i.e. the macrobases) on scalability. First, we consider a 3D example in order to help show that scalability is verified when the micro time scale and the macro time scales are identical. However, when the temporal macrobasis defined in [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF] is used, scalability is not always completely guaranteed. We present a reasoning which shows that classical (h or p) refinement methods of the temporal macrospace do not solve that issue. Then, we propose a general technique for the enrichment of the temporal macrospace which grants the method robustness.

The spatial macrobasis

No condition is required a priori on the spatial macrobasis defining F M E (see [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF]). However, in order to ensure numerical scalability, the resultant forces and moments at interface Φ EE ′ must belong to the macro part of the forces (see [START_REF] Ladevèze | A micro/macro approch for parallel computing of heterogeneous structures[END_REF]). With such a choice, thanks to Saint-Venant's Principle, the micro complements have only local effect. Therefore, there exists a small macrobasis capable of uncoupling the global and local contributions of each quantity. The resolution of the macroproblem consists in enforcing continuity of the macroforces at each iteration, which propagates the most relevant information throughout the structure. In practice, the macro part of an interface quantity consists of its linear part alone. In 3D, this leads us to consider a macrobasis of dimension 12. Naturally, in the case of a plane interface (and, more generally, of slightly curved interfaces), the three out-of-plane modes can be disregarded, which reduces the dimension of the basis to 9 (see Figure 1).

In order to illustrate the fact that with respect to space such a macro separation is optimal, as proved in [START_REF] Ladevèze | A micro/macro approch for parallel computing of heterogeneous structures[END_REF], let us consider the evolution over [0, 10s] of the heterogeneous cube-shaped structure of Figure 2, made of two Maxwell-type viscoelastic materials with Young's moduli E 1 i, Poisson's ratios ν i and viscosities η i . The constitutive relations are such that B i = 1 η i K i . The structure is made of a matrix (E 1 = 50GPa, ν 1 = 0.3 and η 1 = 10s) with fiber inclusions (E 2 = 250GPa, ν 2 = 0.2 and η 2 = 1000s). This structure, subjected to prescribed forces F d on its upper face and clamped at its base U d (Figure 2), was meshed with hexahedra, leading to a problem of approximately 100, 000 DOFs. Three studies were performed, in which the structure was divided into 8, 27 and 216 subdomains respectively (see Figure 3).

In order to emphasize the optimality of the spatial basis presented before, we first studied a time-dependent problem using a computational strategy which was multiscale only in space. In order to do that, we divided the time interval [0, 10s] into 20 coarse intervals, each of which consisted of a single micro subinterval I C i = I f j . The temporal macrobasis was defined as a constant function over I C i , so the macrodescription and the microdescription would match. In the following development, this multi-space-scale, single-timescale approach will be denoted SMu.

The problem of Figure 2 was solved with SMu using the three different decompositions of Figure 3. The corresponding evolutions of the convergence indicator η presented before are shown in Figure 4. One can observe that the curves are identical, which means that the convergence of SMu is independent of the decomposition and, therefore, is scalable. The choice of the spatial macrofields has a clear mechanical meaning due to the fact that Saint-Venant's principle ensures the uncoupling of the different scales with a reduced basis.

Remark. When the number of subdomains increases, the size of the macroproblem, which is global with respect to the entire set of interfaces, also increases. Scalability could be affected if the resolution cost of the macroproblem became prohibitive. For such situations, alternative techniques for the resolution of the macroproblem involving the introduction of a third scale have been developed [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF][START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminate[END_REF].

Choice of the temporal macrobasis

When using the multiscale-in-space-only (SMu) approach, each coarse interval consists of a single micro subinterval I C i = I f j and one needs to solve one macroproblem for each interval of the refined time discretization. For time-dependent problems, this stage can represent a significant computation cost, especially if there are many subdomains.

The principle of the time-space multiscale approach, hereafter denoted TSMu, consists in approximating the evolution of the spatial macroquantities over [0, T ]. Then, the macroproblem is no longer solved for each micro time step, but solved "on average" over the coarse grid.

Classically, for the temporal macrobasis, one uses polynomial functions. In order for the macroproblem to remain small, that basis is usually limited to quadratic functions [START_REF] Ladevèze | On a multiscale computational strategy with time and space homogenization for structural mechanics[END_REF] (see Figure 5). In practice, the temporal macrobasis is orthonormalized in the sense of the scalar product:

a(t), b(t) I C i = I C i a(t)b(t)dt (9)
The problem of Figure 2 was solved using TSMu. In order to do that, we considered that the time tical seems to indicate that in this case TSMu is scalable. The major advantage of this approach is that at each iteration the resolution of 20 global problems homogenized only in space is replaced by the resolution of a single global time-space homogenized problem without affecting the convergence of the method. According to the tests we performed, this scalability property is often verified, but a partial loss of scalability has been observed in some more complex cases. The objective of this paper is to pinpoint this problem and attempt to provide a means to circumvent it.

f 1 M t i C t i+1 C f 2 M t i C t i+1 C f 3 M t i C t i+1 C t t t

The partial scalability issue

The part of the forces which, according to Saint-Venant's Principle, propagate throughout the whole structure is associated at each time step with the resultant and moment of the forces at the interface. Due to the fact that the macro part of the forces is an average in time, the resultant and moment of the micro part are nonzero at the interfaces and need to propagate throughout the structure. This cannot be taken into account in the homogenized problem and must be dealt with through a single-scale treatment, which results in a partial loss of scalability. Most of the time, when the macrospace is well-adapted to the problem, as in the previous example, this phenomenon does not occur.

Here, in order to highlight the problem using a relatively artificial problem in which the time evolution of the prescribed force is not represented exactly by the temporal macrobasis of TSMu, let us consider the evolution of the viscoelastic heterogeneous structure of Figure 7(a). The structure is clamped at the base and subjected to a traction force at the top. The materials are the same as in the previous example. The time interval is divided into 60 microintervals I f j . Figure The problem was solved using successively SMu, TSMu and the full single-scale approach, denoted Mo, in which the method is applied without solving a homogenized problem. Since SMu has been proven to be optimal and scalable [START_REF] Ladevèze | A micro/macro approch for parallel computing of heterogeneous structures[END_REF][START_REF] Ladevèze | A micromacro and parallel computational strategy for highly heterogeneous structures[END_REF], and since the homogenized problem of TSMu is an approximation of the problem of SMu, SMu will be used as the reference for TSMu.

Figure 8 shows the convergence curves of the three methods. One can observe that even though TSMu shows significant improvement compared to the single-scale approach Mo, it suffers from a significant loss of scalability, as quantified by the distance between its convergence curve and that of our reference SMu. Indeed, the oscillations which are part of the loading are not taken into account by the macroproblem and propagate thanks to the single-scale process alone, which is known not to be scalable. This artificial example shows that in some cases TSMu can become only partially multiscale, which is a drawback regarding the robustness of the method.

Later, we will attempt to develop a technique for making TSMu scalable, but first let us study some classical refinement techniques.

Classical refinement techniques

In this section, we will study some simple refinement techniques which consist in enriching the space of the temporal macrofields while keeping the target microscale fixed. Therefore, the converged solution remains the same, and only the convergence rate changes. There are two possible options: a p-refinement option, which consists in increasing the degree of the polynomial functions of the temporal macrobasis (i.e. increasing the frequency bandwidth that the temporal basis can take into account); and an h-refinement option, which consists in decreasing the size of the support of these functions (i.e. increasing the number of coarse intervals I C i making up I = [0, T ]).

The p-refinement option

Here we will consider only the case of polynomial functions, but it would be possible to use more complex families of functions, such as Fourier series or wavelet functions. Problem 7 was solved using TSMu with different polynomial degrees p of the temporal macrobasis, starting from the previous quadratic basis p = 2 and increasing the degree p gradually (p = 3, 5, 7, 9, 11 and 13) until the method became close to SMu. Figure 9 shows the corresponding convergence curves (in black) superimposed on those of Figure 8 (in grey). One can see from that figure that a polynomial basis of at least degree p = 13 is required to come close to the convergence of the reference method SMu. In 3D, the use of such a technique would lead for each interface to a basis of dimension 9 (for the spatial functions) ×14 (for the temporal functions) = 126, compared to 27, which would lead to a much too expensive macroproblem. That, plus the fact that for a more complex problem degree 13 would certainly be insufficient, makes that option clearly inappropriate.

The h-refinement option

Now, let us set the degree of the polynomial functions equal to p = 2 and look at the effect of a refinement of the coarse grid. We divided the time interval [0, T ] into n h = 1 coarse intervals I C i , then gradually increased n h (n h = 2, 3, 6 and 12) until the method came close to SMu, keeping the total number of microintervals I f j the same. error indicator η

n h =1 n h =2 n h =3 n h =6 n h =12
Figure 10: h-refinement applied to the temporal macro approximation requires a rather large number of coarse intervals n h = 12, and the associated cost is close to that of SMu. In this case, the microdescription and the macrodescription are quasi identical. However, contrary to the p-refinement approach, the advantage of this technique is that it converges because the asymptotic case is quasi equivalent to the reference method SMu.

In conclusion, these simple refinement techniques are unsuitable because they require such a level of enrichment that the computation costs become prohibitive. In addition, these techniques are absolutely not predictive because parameters h and p must be known a priori and there is no criterion which can help with their choice.

In the next section, we propose an automatic and adaptive enrichment procedure which overcomes that difficulty.

An automatic enrichment technique

In this section, we describe a technique which consists in enriching the temporal macrobasis automatically in order to make it as compatible as possible with the linear part (i.e. the spatial macro part) of the interface's force and velocity fields for each substructure.

Extraction of the enrichment functions

At the end of the linear stage of an iteration, the projections of the microresiduals F m EE (x, t) and Ẇ m EE (x, t) onto the spatial macrobasis functions e M k for each microinterval I f j are extracted. These projections, denoted r F (t) and r W (t), are defined by:

r F (t) = Φ EE ′ e M k (x) • F m EE ′ (x, t)dΦ r W (t) = Φ EE ′ e M k (x) • W m EE ′ (x, t)dΦ
These time functions, on the refined temporal grid, correspond to the evolution of microphenomena which have nonzero resultants and moments (or rigid body modes) at the interface. According to Saint-Venant's principle, these microphenomena have a global effect, yet they are not taken into account by the macroproblem. Thus, the macro procedure is not applied to the entire resultant and moment forces at the interfaces, which is the reason why scalability is partially lost.

The principle of the method proposed here is to enrich the temporal macrobasis with these residual functions. However, these functions are not sufficient to restore scalability. An important condition is also that the enrichment of the temporal macrobasis be consistent with the material model being used. The temporal macrobasis is common to both forces and velocities. If a function r F (t) extracted from the force fields is capable of enriching the macrobasis and improving the representation of the evolution of the resultant forces and moments in the solution, the same function is not necessarily capable of improving the representation of the evolution of the corresponding displacements. For example:

• For a linear elastic material, stress and strain are proportional:

σ = Kε
Then, a function r F (t) is sufficient to improve the representation of both stress and strain.

• for a viscoelastic material such as in our examples, the model is based on a duality between the strain rate and a combination of the stress and its time derivative:

ε = K -1 σ + B(σ)
Thus, when a time function capable of enriching the temporal macrobasis is extracted from the stresses, its time derivative is added, too. More generally, the derivatives of the functions in the temporal macrobasis must also belong to that macrobasis. This makes the initial choice of a basis of polynomials a good choice.

Remark. Here, for the sake of simplicity, we chose to use the same temporal macrobasis for the velocities and for the forces. One could consider using different bases. This would certainly lead to a reexamination of our choice: for example, one could enrich the force basis using one function and the velocity basis using a sort of "dual" function. This possibility needs to be studied in a future work.

Selection of the most relevant functions

If n M esp denotes the number of spatial macrofunctions and n Φ the number of interfaces, there are m = 2 × n M esp × n Φ such functions which can be added to the temporal macrobasis. Of course, should the temporal macrobasis be enriched using all these residual functions, the size of the macroproblem would explode and its resolution cost would be prohibitive. Therefore, it is necessary to find a way to rank these functions in order to select only the most relevant ones.

The objective is to find a relatively small number of functions whose linear combination could approximate the whole set of m functions as well as possible. A convenient way to do that consists in calculating a truncated Proper Orthogonal Decomposition 2 (POD, see [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF]) of these functions.

Let n f i be the number of refined subintervals of the coarse interval I C i . We now have at our disposal m residual functions {r j (t)} 1 j m (r in their discrete form) arranged in a m × n f i matrix A, such that A ij is the value of the j th function for the i th refined subinterval.

A = . . . r j . . .
The calculation of a truncated proper orthogonal decomposition of A comes down to the calculation of the first few eigenfunctions of the covariance matrix C of A (see [START_REF] Chatterjee | An introduction to the proper orthogonal decomposition[END_REF]):

C = AA T
Let us note that in our case, since n f i ≪ m, C is a small n f i × n f i symmetric matrix. The calculation of the eigenvectors of that matrix is relatively inexpensive and independent of the number of residual functions m. Let us sort the eigenvalues λ i in ascending order:

|λ 0 | . . . |λ i | . . . |λ n f i |
The eigenfunctions associated with the first n eigenvalues (λ i ) 0 i n constitute the best basis of size m, i.e. one whose linear combination leads to the best approximation of the complete m functions in the sense of the Frobenius norm, which is a discrete version of the L 2 norm:

A 2 f = tr(AA T )
2 Also known as Singular Value Decomposition, Principal Component Analysis or Karhunen Loeve Expansion, see [START_REF] Liang | Proper orthogonal decomposition and its applications-part I: Theory[END_REF] This means that the proper orthogonal modes are optimal and orthogonal in the sense of the L 2 norm. Let us recall that the macrobasis must be orthonormal in the sense of the scalar product •, • I C i defined by Equation ( 9), which is more relevant physically. Therefore, the proper orthogonal modes should be optimal in the sense of the following metrics:

r i (t) 2 I C i ,m = m i=1 r i (t), r i (t) I C i (10)
In the case of heterogeneous partitions of I C i , the calculation of the POD must be modified slightly. Let us introduce the n f i × n f i diagonal matrix D whose elements are equal to the sizes D jj = t f j+1t f j of the refined intervals I f j . The POD must be carried out in the sense of the following modified Frobenius norm:

A 2 mf = tr(A T DA) = tr m i=1 r T i Dr i = m i=1 tr(r T i Dr i ) = m i=1 r T i Dr i
which is the the discrete version of Norm [START_REF] Combescure | A numerical scheme to couple subdomains with different timesteps for predominantly linear transient analysis[END_REF].

Calculating the POD with such a norm is equivalent to seeking the eigenvalues and eigenvectors of the following modified covariance matrix C: C = AA T D Indeed, let us consider a first-order POD approximation λuv T of A, where λ is a scalar and u ∈ R n f i and v ∈ R m are normalized vectors. λ, u and v are unknowns, which are sought in order to minimize the following metrics:

A -λuv T 2 mf = tr(A T DA) -2λtr(A T Duv T ) + λ 2 tr(vu T Duv T )
We seek v which minimizes the previous metrics over R m , leading to ∀v ⋆ ∈ R m :

-λtr(A T Duv T ⋆ ) + λ 2 tr(vu T Duv T ⋆ ) = 0 which is equivalent to:

A T Du u T Du = λv (11)
We proceed in the same way with u, which leads to:

Av v T v = λu (12)
Equations ( 11) and ( 12) yield:

AA T Du v T v u T Du = λ 2 u ( 13
)
Let us recall that u and v are normalized, i.e., for example, v T v = u T Du = 1. Equation ( 13) is an eigenvalue problem of C = AA T D in which u is an eigenvector and λ 2 the corresponding eigenvalue. Finally, v can be calculated a posteriori thanks to [START_REF] Cresta | Nonlinear localization strategies for domain decomposition methods in structural mechanics[END_REF]. A study of the properties of such an eigenvalue problem is given in [START_REF] Ladevèze | Nonlinear computationnal structural mechanics-New approaches and non-incremental methods of calculation[END_REF].

The resulting eigenfunctions {u i } 1...n f i of C are orthonormal in the sense of:

r i , r j = r T i Dr j (14)
Indeed, let (u 1 , u 2 ) ∈ (R n f i ) 2 be two eigenvectors of C associated with the eigenvalues (λ 2 1 , λ 2 2 ) ∈ R 2 such that λ 1 = λ 2 . Then, one has:

u T 1 DAA T Du 2 = λ 2 1 u T 1 Du 2 = λ 2 2 u T 1 Du 2 Since λ 1 = λ 2 , one has: u T 1 Du 2 = 0 (14)
is the discrete version of the scalar product [START_REF] Cognard | A large time increment approach for cyclic plasticity[END_REF] which enables one to add the eigenfunctions u directly to the temporal macrobasis.

As we will see in the example, the eigenvalues of C (or C) decrease rapidly, so the first few eigenfunctions are the dominant ones. For example, one can use the following indicator to select the most relevant functions:

|λ i | |λ 0 | < ǫ λ (15)
These eigenfunctions are suitable functions to add to the temporal macrobasis. In order to keep the macroproblem small, one must limit the enrichment to the first few functions, but thanks to the POD, these first functions are the most relevant ones.

Back to the example

This automatic enrichment technique (which will be denoted aTSMu) was applied to the previous problem of Figure 7. The temporal macrospace was characterized by a single macrointerval and the initial basis consisted of only 3 quadratic polynomials. The number of functions which could be added to that basis was limited to 3. This stage was carried out every 10 iterations and the truncation criterion was chosen to be equal to ǫ λ = 10 -4 . The corresponding convergence curve (in black) is shown in Figure 11. One can observe on that figure that with only a few additional functions the convergence can be accelerated significantly. These additional functions, which do not require any a priori information and are calculated automatically during the iterations, grant the method a convergence rate which is close to that of SMu, which is scalable.

Remark. Since the temporal macrobasis has been modified, the homogenized operators need to be updated. The associated computation cost is significant, but one does not need to use this automatic enrichment procedure at each iteration. Indeed, in Figure 11, the discontinuous black curve corresponds to one enrichment per iteration, while the continuous black curve corresponds to one enrichment every ten iterations. The fact that these two curves are similar means that the frequency with which the basis is updated influences the effectiveness of the technique only slightly. Thus, the cost of the process can be easily reduced. Furthermore, in the case of nonlinear behavior, the homogenized operators need to be updated from time to time as well because of the modification of the search directions. In that case, one should update the temporal macrobasis and the search directions simultaneously, which makes the additional cost of the automatic enrichment technique very affordable.

When the macrospace is compatible with the data

The objective of the previous example was to make the procedure easily understandable, but the problem is somewhat artificial. In general, it is necessary to verify that the temporal macrospace (the macrobasis and the h and p parameters) is compatible with the data. In particular, the entire resultant forces and moments of the loading must belong to the macro part. In the previous example, this would have consisted in introducing, for example, the sinusoidal functions of the perturbation into the macrobasis. Even in this case, the evolution of the resultant forces and moments at the interfaces may be very different from that of the loading, for example in the context of nonlinearities.

Let us consider the problem of Figure 12(a) where the same structure as before is subjected to a lateral force (left) in the presence of a rigid solid (right). The rigid solid was modeled using an interface with special frictionless contact behavior and a gap g = 1% of the size of the domain. Since the time evolution of the loading (Figure 12(b)) is quite simple, the quadratic temporal basis was compatible with the data. Figure 13 shows the convergence of methods SMu, TSMu, Mo and aTSMu using the same parameters as in the previous example.

However, although the macrospace was compatible with the data, it was insufficient to take into account the temporal phenomena in the vicinity of the contact. These phenomena were dealt with in the single-scale process, which resulted in a degradation of the overall convergence of the method. Nevertheless, in that example, thanks to the automatic enrichment phase, aTSMu converged very rapidly and required 10 times fewer iterations than the standard TSMu to bring the error down to η = 10 -3 , which corresponds to a convergence rate comparable to that of the reference SMu.

A three-dimensional example with frictional contact

In order to illustrate the strategy on a more difficult example, let us consider the quasi-static evolution over [0, 10s] of the viscoelastic 3D structure of Figure 14 containing a crack over 4/5 th of its length. The structure is made of a Maxwell- type viscoelastic material whose constitutive relation is such that B = 1 η K, where K is Hooke's tensor with Young's modulus E = 100GP a, Poisson's ratio ν = 0.3 and viscosity η = 10. The structure is clamped at the base and subjected to a preloading force at the top and a traction force along the upper lip of the crack. The domain was meshed with hexahedra with approximately 100, 000 DOFs. The time interval consisted of a single coarse interval divided into 60 refined subintervals. The domain was divided into 10 substructures and 47 interfaces. Using a classical choice for the macrobases (linear in space and quadratic in time) the macrospace was compatible with the data. The crack was modeled by a frictional contact interface whose formulation can be found in [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF].

The loading was chosen such that all possible scenarios would be encountered at the contact interface (see Figure 15 was solved using Mo, SMu, TSMu and aTSMu successively. In the case of aTSMu, the enrichment phase was carried out every 10 iterations, the truncation criterion was taken to be equal to ǫ λ = 10 -4 and the maximum number of additional functions in the basis was set to 3. Figure 16 shows the corresponding convergence curves.

In this example, one can observe that the stan- dard time-space multiscale approach TSMu was penalized by the time approximation of the evolution of the resultant forces and moments to the point that the convergence indicator was similar to that of the full single-scale approach Mo. As in the previous example, the automatic enrichment technique made the convergence of aTSMu similar to that of the reference SMu, in which the temporal macrospace contains the entire fine description and which is scalable. Thus, the proposed approach seems to be robust even in the case of more complex problems. It also seems to give the time-space approach quasi full scalability thanks to the adaptation of the temporal macrobasis.

Conclusions

In this paper, using several examples, we showed that in some cases the seemingly most natural separation a priori of the temporal scales-slow scale / rapid scale -is not always optimal and leads to a nonscalable LATIN algorithm. Refinement techniques of the h-type (which refine the temporal macrospace) and of the p-type (which introduce higher frequencies) were found to be relatively ineffective and, moreover, too problemdependent. Therefore, we raised the question of the separability of the scales. We proposed a vision which consists in considering the macro part of a quantity to be its projection onto a reduced basis. This choice leads to a smaller temporal macrospace, which is well-adapted to the problem and its loading. This reduced basis, which is unknown a priori, is adapted automatically during the iterations of the calculation and leads to a negligible increase in resolution cost in the general case of nonlinear evolution laws.

The technique consists in selecting a set of temporal residuals which reflect the high-frequency phenomena which have a global influence on the spatial level. Since these phenomena are microscopic, they cannot be taken into account in the macroproblem, which affects the scalability of the method. However, a technique based on the POD approach enables one to select only the dominant functions from that set, leading to an enrichment of reasonable size. This time-adaptive separation leads to a numerically scalable calculation method regardless of the frequency content of the loading.

One of the consequences of this work could be the modification of the current algorithm, leading to a macroproblem which could be solved using Proper Generalized Decomposition [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF][START_REF] Chinesta | Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization[END_REF][START_REF] Nouy | A generalized spectral decomposition technique to solve a class of linear stochastic partial differential equations[END_REF]. Thus, the reduced basis could be generated a priori at each iteration without resorting to POD. This would also enable the extension of this type of reduced-basis micro/macro decomposition to the spatial variable, which could probably lead to a reduction of the spatial macrospace, which in some three-dimensional cases can be more extensive than strictly necessary.

→

  The macroproblem leads to ˙ W M and, through another set of micro resolutions (3), to s 2 E . Then, one can determine s E completely. Initialization: s 0 ∈ A d while η tol. do • Local stage: ŝ ∈ Γ Nonlinear problems solved locally at each point of the discretization of Ω and I. • Possible recalculation of the search directions • The linear stage: s ∈ A d → The microproblem (1st set) linear predictions defined independently over Ω E × I C i , ∀E. → The macroproblem a time-space homogenized problem defined over Ω × I C i The micro problem (2nd set) linear predictions defined independently over

Figure 1 :

 1 Figure 1: The spatial macrobasis

Figure 2 :

 2 Figure 2: Definition of the problem

Figure 3 :

 3 Figure 3: the different decompositions of the structure

Figure 4 :

 4 Figure 4: Convergence of SMu for different domain decompositions

Figure 5 :

 5 Figure 5: The current choice of the temporal macrobasis

Figure 6 :

 6 Figure 6: Convergence of TSMu and SMu for different domain decompositions

Figure 7 :

 7 Figure 7: A problem with a macrobasis incompatible with the loading 7(b) shows the evolution of the prescribed force

F

  d and its projection F M d onto the temporal macrobasis of TSMu. One can immediately see that the resultant of the micro part of the loading is nonzero, which leads one to imagine that the homogenized problem of TSMu cannot propagate the loading properly.

Figure 8 :

 8 Figure 8: Convergence of the SMu, TSMu and Mo methods

Figure 9 :

 9 Figure 9: p-refinement applied to the temporal macro approximation

  Figure 10 shows the corresponding convergence curves (in black). As with the previous case, this technique

Figure 11 :

 11 Figure 11: Convergence of the aTSMu method

Figure 12 :Figure 13 :

 1213 Figure 12: The frictionless contact problem (with a macrobasis compatible with the data)

Figure 14 :

 14 Figure 14: Definition of the problem

Figure 15 :

 15 Figure 15: The solution of the frictional contact problem

Figure 16 :

 16 Figure 16: Convergence curves for SMu, TSMu, aTSMu and Mo

The notation E denotes the restriction of a quantity to substructure Ω E