A separated representation involving multiple time scales within the Proper Generalized Decomposition framework

Solutions of partial differential equations can exhibit multiple time scales. Standard discretization techniques are constrained to capture the finest scale to accurately predict the response of the system. In this paper, we provide an alternative route to circumvent prohibitive meshes arising from the necessity of capturing fine-scale behaviors. The proposed methodology is based on a time-separated representation within the standard Proper Generalized Decomposition, where the time coordinate is transformed into a multi-dimensional time through new separated coordinates, each representing one scale, while continuity is ensured in the scale coupling. For instance, when considering two different time scales, the governing Partial Differential Equation is commuted into a nonlinear system that iterates between the so-called microtime and macrotime, so that the time coordinate can be viewed as a 2D time. The macroscale effects are taken into account by means of a finite element-based macro-discretization, whereas the microscale effects are handled with unidimensional parent spaces that are replicated throughout the time domain. The resulting separated representation allows us a very fine time discretization without impacting the computational efficiency. The proposed formulation is explored and numerically verified on thermal and elastodynamic problems.


Introduction
Many engineering problems are defined in very large time intervals (e.g., when dealing with fatigue, aging, dynamics with loadings involving multiple characteristic times) and, at the same time, the response must encompass the different time scales present in the model. Following standard time marching approaches, a suitable time step that captures the evolution of the finest scale has to be adopted to ensure a reliable modeling, deriving into a prohibitive simulation cost. Many numerical algorithms have been proposed to bypass the time marching approach up to the finest time scale.
For instance, when only macro-scale effects are envisaged, time-homogenization techniques perform well, although subject to a number of hypotheses which are not always realistic, as well as to the scale separation.
In [1], it is proposed a technique that makes use of the separation of variables. The time domain is partitioned into a number of subdomains, where times defining the subdomain partition describe the macrotime, whereas the subdomain discretization resolves fast responses (microtime). The main issue characterizing this technique is the necessity of ensuring the continuity of the resulting two-scale discretization. Lagrange multipliers are adopted for this purpose. This choice makes substantially prohibitive the computational implementation associated with this procedure.
In [2][3][4] a numerical algorithm based on the LArge Time INcrement-Proper Generalized Decomposition (LATIN-PGD) is presented. Temporal and spatial multi-scale behaviors in solid mechanic problems are handled by associating a temporal macro basis with the interfacial degrees of freedom coupling different macro domains. The scalability of the methodology is restored via an appropriate correction of the temporal basis founded on the residual.
Also in [5], the time domain is partitioned into subintervals, and a common reduced basis is adopted in every subinterval. Afterwards, additional interface restrictions are imposed to ensure the continuity of the primal variable and of the associated time derivative. Hence, values at the endpoints of each macro interval are enforced to vanish, and the offset is calculated on the macro discretization.
An alternative view is proposed in [6,7], within the so-called Parareal scheme, that iterates between macro and micro domain partitions, so that the initial conditions for each micro interval are provided by the macro resolution of the problem. The algorithm is easily parallelizable, making it very efficient for addressing general models.
Another numerical strategy based on combining different PGDs under the partition of unity (PU) rationale is proposed in [8,9]. Here, the main idea is to use macro shape functions fulfilling the partition of unity. This ensures moving throughout different PGDs in diverse intervals while guaranteeing a perfect continuity. Actually, the overlap between the PGDs solutions related to the overlapping subdomains ensures the continuity when the PGDs solutions are multiplied by the macroscopic shape functions thanks to the standard PU features. Despite its effectiveness, the necessity of combining microscopic discretization with the macroscopic PU enrichment, hinders the computational implementation of this approach.
This paper aims at computing the time evolution of the unknown fields involved by a multiscale Partial Differential Equation (PDE), taking into account all the information that it involves, i.e., covering all the time scales. The solution procedure that we are proposing includes neither major hypotheses, nor the need of time scale separation, as for time homogenization techniques. In more detail, the new approach ensures continuity in a direct and computationally inexpensive manner without resorting to the use of Lagrange multipliers, penalty or the PU paradigm. Our proposal directly derives from the discrete tensor formulation of the separated representation involved in the PGD constructor, making use of a tensor formalism. For that purpose, in the next section, after providing the basics for a PGD tensor representation, 1 we introduce the new approach to manage temporal multi-scale problems and we apply it to the heat equation. In "Multiscale resolution of transient elastodynamic problems" section we extend such an approach to computational solid mechanics and, more concretely, to elastodynamics. Finally "Conclusion" section shortly summarizes the main achievements of the present work.

Tensor-based PGD methods
Here we shortly recall the PGD algorithm in a general tensor framework, as described in [13,14]. For a detailed formalization of the PGD method, we refer the interested reader, e.g., to [15][16][17][18][19][20]. The solution of a multidimensional problem is built as a sum of tensor products of functions defined in some sub-spaces with moderate dimension (1, 2 or 3), thus providing a general separable representation form. In particular, the authors make use of the best rank-1 approximation property of tensors of order 3 or higher [21] (in [22], it has been proved that tensors of order 3 or higher can fail to have best rank-n approximation) to propose an iterative method based on the so-called projection-enrichment technique [13].
Let be a multidimensional domain involving several coordinates x i (not necessarily one-dimensional), which can coincide, for instance, with spatial coordinates, time, model or geometric parameters. We consider the weak formulation of a linear problem: with V ( ) an adequate function space ensuring the well-posedness of such a formulation. We assume that, after a discretization of problem (1), we are lead to solve the linear system: where the operator on the left-hand side and the right-hand side member are expressed in a separated form as The approximated PGD solution ψ is sought in the discrete separated form and an analogous representation is adopted for the test function ψ * , being Solution ψ is built-up by using a projection-enrichment iterative scheme [13]. In particular, it will be assumed that the global convergence is attained when the error estimator ε = ||A ψ − B|| 2 is small enough, where ||·|| 2 denotes the standard Frobenius norm of a tensor of order d.
The projection stage consists of finding the set of coefficients α j in (4) verifying the relations where The enrichment stage includes new candidates for enriching the reduced separated approximation basis, so that ψ can be updated as and (2) is replaced by the system In particular, within a fixed point alternating direction algorithm, at each iteration we look for the computation of a single discrete function, r j , all the other components of ψ r being assumed known, and after setting ψ * to r 1 ⊗ · · · r j−1 ⊗ r * j ⊗ r j+1 · · · ⊗ r d . This strategy leads to solve the linear system where For the explicit computations leading to system (8), we refer the interested reader to [13], whereas a convergence analysis for this greedy rank-1 update algorithm can be found in [21].
Remark 0.1 The projection in (5) does not represent the only possible choice when defining ψ. Alternatively, one could use a standard greedy-based enrichment [14,21]. The aim of the projection is to regularize the solution, thus getting rid of possible spurious modes yielded by an enrichment procedure and, consequently, to enhance the global convergence. Indeed, since modes are normalized, coefficients α j in (5) can be intended as a truncation error since they express the importance of the different enrichment terms. This allows us to stop the enrichment procedure when not significant modes occur, thus limiting the inclusion of noisy information.
Remark 0.2 As discussed in [13], the convergence of the fixed point strategy characterizing the enrichment step is guaranteed for symmetric discrete operators. Actually, numerical tests performed with non symmetric discrete operators (e.g., with hyperbolic operators) exhibit some difficulties in the convergence of the fixed point procedure as well as a loss of optimality in terms of number of function products to be computed. For this reason, all the problems analyzed in this work are suitably rewritten in order to deal with a symmetric operator. In general, the new problem to be solved reads: and the associated projection stage becomes Analogously, system (7) will be replaced by the new one so that system (8) takes the form

A two-time scale separated representation within the PGD method
In the proposal of numerical methods to manage multi-scale phenomena, a crucial point is represented by the capability to ensure the continuity of the approximation. This turns out to be an issue also when dealing with a PGD approach. For instance, in [9], PGD is combined with the partition of unity (PU) paradigm to circumvent the continuity problem.
Thus, a generic multiscale function u is expressed as where the collection, , of the so-called centroids x i introduces a macro-partition of the time window, N i (x) is the macroscale shape function associated with the i-th centroid, x i , ξ (x − x i ) is a dependent variable characterized by an offset based on x i , and G j (ξ ) denotes the j-th microscale shape function, with g j the corresponding microscale degree of freedom. Here, the microscale effects occurring in the compact support of the macroscale shape function N i (x) are mapped into the corresponding microscale space, via ξ . Formulation (9) allows the authors to automatically ensure the continuity to the approximation u as proved in [9]. Other formulations enforce the continuity in different way, for instance by resorting to Lagrange multipliers as in [1].
The present work takes a step forward with respect to the available literature by proposing a new approach which ensures the continuity of the PGD solution in a computationally inexpensive manner, without resorting to Lagrange multipliers or penalization strategies.
We exemplify the new method on a simple problem, i.e., on the heat equation solved in a one-dimensional spatial domain x = (0, L) over the time window t = (0, T f ), and completed with non homogeneous Dirichlet data, with k the diffusivity coefficient and f = f (x, t) the source term.
A standard time-marching approach requires a discretization of t by means of N t times. On the contrary, the technique here presented is based on splitting the whole time i=0 of equispaced macro-dofs and by setting T 0 = 0 and T M = T f . Afterwards, in the first macro subinterval, a new continuous time-microscale, τ ( 1 ), is introduced so that the subinterval 1 is discretized by involving m + 1 micro-dofs, {τ 1 j } m j=0 and setting τ 1 0 = T 0 . Also the discretization at the microscale is uniform, with a constant time step set to τ throughout 1 . Moreover, the last microdof is internal to 1 , being τ 1 m = T 1 − τ . By means of a tensor product, the micro discretization in 1 is replicated throughout the other subintervals i+1 , covering the whole time window (see Fig. 1 for a sketch). In such a way, a fine time grid consisting of N t = M(m + 1) dofs is obtained by mapping m + 1 micro-dofs onto a coarser grid of M macro intervals, without a priori needing a fine discretization of the whole time window. In particular, the value of the solution at the right endpoint, T i+1 , of each subinterval i+1 is recovered by means of this tensor product. As a consequence, such a discretization is intentionally not meant to predict values at T M = T f . Figure 2 shows the Cartesian product of the macro and micro dofs together, with the associated nodal evaluations of the generic unknown function u in (10).  If we discretize the time derivative in (10) by means of a standard finite difference formula, we have: that we can write also as with where the first equation understands that the initial value u 0 has to be included in the right-hand side of the equation associated with n = 0, while the subscript N t for the operator E denotes the size of the corresponding matrix. This notation will be particularly useful to simplify the discussion which follows. Using a two-scale time variable, the time derivative in (11) is replaced by the new one i=0 is adopted to identify the subinterval i+1 of interest, while τ ∈ [T 0 , T 1 ) represents the actual temporal independent variable, spanning 1 , i.e, via tensor product, the generic subinterval i+1 . This leads to perform the time derivation with respect to variable τ .
With a view to the application of PGD for approximating problem (10) with the twoscale time discretization, matrix E N t in (13) associated with the standard time coordinate deserves to be modified in order to replicate the PGD tensor structure. The idea is to introduce a suitable macro-micro decomposition of the time derivative operator in (12). For this purpose, let us identify the partial derivative of order zero with respect to T , ∂ 0 ∂T 0 , with the identity operator, so that, ∂ 0 ∂T 0 f (T ) = f (T ), and, as a consequence, identity holds true. This equality suggested the authors in [1] to propose the approximation for the standard time derivative, which will be investigated and further developed throughout this work. Approximation (14) can thus be written in matrix form as where now the subscript M refers to the number of blocks constituting matrix E T 0 ,τ , while matrix E m+1 , sized according to discretization at the micro-scale, is defined as in (13), after replacing t with τ , and replicated throughout the subintervals i+1 by means of the Kronecker product with the discrete identity operator I M . 2 2 We observe that symbol ⊗ refers to a revision of the standard Kronecker product, being meant However, definition (15) leads to have some entries along the first sub-diagonal which are identically null. This happens essentially when switching from a block to the subsequent one, namely, when moving from the last internal dof associated with the micro-scale in a certain subinterval, i+1 , to the first dof associated with the macro-scale of the subinterval i+2 , resulting in a loss of continuity between the fine and the coarse time grids. An effective, yet computationally demanding, remedy to this issue is based on enforcing the continuity in all the macrodofs, either by using a penalized formulation or by employing Lagrange multipliers, as detailed in [1].
In this work, we propose to exploit the properties of the tensor product to recover the equality in (15). To this aim, we introduce the decomposition which, in a compact form, reads as where matrices C m+1 and L M are defined as respectively with δ ij the Kronecker symbol. The decomposition in (17) allows us to reconstruct the correct time derivative operator while preserving the macro-micro decomposition characterizing the two-scale time derivative operator. This method embraces the PGD formalism while reducing significantly the computational cost, since no extra computational effort is required to enforce continuity at each dof associated with the macroscale.
Remark 0.3 The multiscale approach proposed in this paper exhibits two main benefits: 1. The possibility of describing in an efficient way physical phenomena with the fine scale resolution; 2. A description richness equal to N t = M(m + 1) accomplished with a computational complexity of the order of M + m + 1, thus yielding considerably reduced memory allocation requirements to store time operators.
Footnote 3 continued respectively. This is a common choice in the PGD community to preserve the same order of the variables in the separated form [23].
Coming back to the approximation of the heat equation in (10), a classical discrete PGD approximation, based on a space-time separated form, reads: where N U denotes the number of the PGD modes, index 1 refers to the space dependence, index 2 is associated with the time, B H is a suitable discrete separated approximation of the source term f (t, x) evaluated at some discrete points, A H is a discrete separated representation of the differential operators in (10), built by exploiting the tensorial identity and assembled according to the discretization adopted for the space (e.g., finite differences, finite elements) and for the time. With the same notation as in (14), we identify the partial derivatives of order zero with the corresponding identity operator. For instance, let us assume to discretize (10) by employing the implicit Euler scheme for the time and a centered finite difference approximation for the space derivative, with N x the number of spatial dofs and N t the number of discrete times. Then, the operator A H can be separated as where E N t denotes the matrix associated with the implicit Euler scheme, I N x and I N t are the identity matrices in space and time, respectively and D N x is the matrix of the discrete Laplacian. Figure 3 shows the sparsity pattern associated with four matrices A k i . Now, in order to apply the new PGD formulation based on a two-scale time discretization to problem (10), we have to recast the differential operators in (20) into a multiscale framework, by adding a further temporal independent variable. This leads us to deal with a framework characterized by two dimensions in time and one dimension in space. For this purpose, we adopt the approximation and we enforce the continuity between the macro and the micro time partitions by resorting to the two-scale decomposition in (17), so that the operator A H in (19) can be separated as where, analogously as in (19), subscript 1 identifies the space dependence, while subscripts 2 and 3 refer to the time microscale and to the time macroscale, respectively. The spatial operators are defined in a straightforward way as whereas the two-scale time discretization leads us to adopt the modified tensorial decomposition in (16), with A 1 2 = E m+1 the matrix characterizing the implicit Euler discretization at the microscale; A 1 3 = A 3 3 = I M the discrete identity matrix associated with the subintervals i+1 ; A 2 2 = C m+1 the correction matrix, sized according to the microscale, with a unique non-null element in the upper-right corner; A 2 3 = L M the lower shift matrix, sized according to the macroscale, with the first subdiagonal of elements all equal to one; A 3 2 = I m+1 the discrete identity matrix associated with the microscale. Figure 4 gathers the sparsity pattern of the nine matrices involved in the separated representation (22).
Concerning the right-hand side in (19), operator B H coincides with a separated representation, with respect to the spatial scale and the two temporal ones, of the source term f (t, x) in (10), evaluated at some discrete points, and computed through a HOSVD [24].

Remark 0.4 A cross comparison between Figs. 3 and 4 highlights that the two-scale PGD
formulation of problem (10) does not affect the spatial discretization. On the contrary, the matrices associated with the two-scale time discretization are much smaller with respect to the ones characterizing the separated representation in (21), with a consequent reduction in terms of memory allocation (in particular, the storage requirement for the time matrices reduces from 3(m + 1)M − 1 to 3(m + 1 + M) − 1, which means of a factor

Multiscale resolution of transient elastodynamic problems
This section generalizes the multiscale PGD formulation detailed on the heat equation in the previous section to a more challenging setting represented by an elastodynamic problem in a multiscale transient regime. In this setting, different strategies based on the separation of the space and time are investigated in [23]. Here, the authors propose a multi-field Proper Generalized Decomposition which adopts standard finite elements to discretize the space, to be combined with several discretization schemes in time, such as the Newmark scheme [25], a single-field time discontinuous Galerkin method, a two-fields time continuous and discontinuous Galerkin approach.
To simplify the discussion, we focus on a transient elastodynamic problem in a onedimensional spatial domain which models, for instance, the traction-compression waves travelling along a linear elastic medium, = (0, L), during time interval I = [0, T ). The scalar displacement field is denoted by u = u(x, t), for x ∈ and t ∈ I. The medium is subject to a displacement, u d = u d (x, t), and to an external force, f d = f d (x, t), applied to portions, ∂ u × I and ∂ f × I, of the boundary of the space-time cylinder, respectively, with ∂ f ∪ ∂ u = ∂ and • ∂ f ∩ • ∂ u = ∅ (notice that, in the specific context at hand, ∂ u and ∂ f reduce to the endpoints of the domain , dealing with a one-dimensional setting). The initial state is described by the displacement field, u 0 = u 0 (x), and by the velocity field, v 0 = v 0 (x), at time t = 0. Finally, the medium is characterized by a density, ρ, an elasticity modulus, E, and a section, A. Thus, the strong formulation of the considered transient elastodynamic problem reads: find u : To simplify the notation, in the following we adopt the conventionu andü to denote the first-and the second-order time derivative, ∂u ∂t and ∂ 2 u ∂t 2 , of u respectively. We employ a standard centered finite difference scheme to approximate the space, after discretizing with N x dofs, so that the associated semi-discrete formulation is: find where M = ρ E I N x and K = − D N x denote the mass and the stiffness matrix, respectively, with I N x and D N x the identity and discrete Laplacian operators, u 0 and v 0 are the vectors sampling the initial data, u 0 and v 0 , at the spatial dofs, and where the boundary data have been properly included in the right-hand side F. To discretize the time, we subdivide the whole time window, I, by means of (N t + 1) uniformly distributed times, t i , such that t 0 = 0, t N t = T , and t i − t i−1 = t, for i = 1, . . . , N t . Following [23], we apply the Newmark scheme to (24), so that we have to solve the fully discretized problem: find u ∈ R N x ×N t , such that where : denotes the two times contracted product between tensors 3 , with 0 ≤ γ ≤ 1 and 0 ≤ β ≤ 1 2 the parameters involved by the Newmark algorithm, a = 1/2 − 2β + γ , b = 1/2 + β − γ , F(t 0 ) is the right-hand side in (24) evaluated at the initial time t 0 , and F collects all the vectors F(t i ) for i = 1, . . . , N t , accounting for the exterior force and the imposed displacement (we refer the interested reader to the reference paper [23] where the explicit computations leading to linear system (25) can be found).
Formulation (25) is instrumental to provide the discrete PGD approximation for problem (23) based on a standard space-time separation, which reads where the same notation as in (19) being adopted for the subscripts identifying the space and the time, respectively, while B E simplifies to after assuming F = 0 (i.e., no external forces) and u 0 = 0 (i.e., a vanishing initial displacement).
The discrete operators involved in the definition of A E have the sparsity pattern shown in Fig. 5, when choosing = (0, 1) and I = [0, 100).
We now generalize the two-scale temporal discretization for a PGD approximation introduced in the previous section to the elastodynamic problem (23). For this purpose, the tensor decomposition correction proposed in (16) for the implicit Euler discrete operator, E N t , is here extended to the Newmark discrete operators, N 1 and N 2 , in (26), so that it holds m + 1 and M denoting the number of the temporal micro-dofs and of the time macro subintervals as in "A two-time scale separated representation within the PGD method" section, with N t = M(m + 1). Decompositions (29) and (30) can be rewritten in a compact form as where matrices N 1,m+1 and N 2,m+1 are defined exactly as N 1 and N 2 , respectively except for the size (equal to m+1 instead of N t ); I M is the discrete identity matrix associated with the macroscale; L M is the lower shift matrix defined as in (18); the correction matrices C 1,m+1 and C 2,m+1 are given by The tensorial decompositions in (31) ensure the continuity between the macro and the micro time scales, analogously to decomposition (17). This allows us to extend the temporal continuity between different scales to A E in (27), after separating this operator in terms of the coordinates x, T and τ , and by exploiting decompositions (31), so that we have subscript 1 referring to the space dependence, while 2 and 3 to the temporal micro-and macroscale, respectively as in the previous section. In particular, the spatial contributions coincide with matrices M and K in (24), being The operators associated with the time microscale are , while the ones related to the temporal macroscale are In particular, the time microscale is characterized by the Toeplitz matrices, A 1 2 and A 3 2 , and by the correction matrices, A 2 2 and A 4 2 , with only three not null elements; the diagonal matrices A 1 3 and A 3 3 and the lower shift matrices A 2 3 and A 4 3 identify the time macroscale. The sparsity pattern associated with the 12 operators in (32) is shown in Fig.  6, for = (0, 1) and I = [0, 100).
The right-hand side in (27) has to be separated, as operator A E , in terms of the spatial scale and of the two temporal ones, so that we have Fig. 6 Sparsity pattern of the discrete PGD operators associated with the separated representation (32), for N x = 51, m + 1 = 400 and M = 5 Remark 0.5 As expected, it turns out that the temporal multiscale PGD discrete formulation does not affect the spatial representation, the sparsity pattern of matrices A 1 1 and A 2 1 in Fig. 5 being preserved by matrices A 1 1 -A 2 1 and A 3 1 -A 4 1 of Fig. 6, respectively. On the contrary, the multiscale time discretization yields matrices with a considerably lower size, albeit in a larger number, with a remarkable saving in terms of memory allocation (notice that the decomposition adopted in the two-scale discretization suits the time grid in

Numerical simulations
In this section we assess in more detail the PGD approximation combining the two-scale Newmark scheme with the centered finite difference method. To this aim, we consider a beam with constant section A = 1, in a traction and compression regime, where an initial velocity, v 0 , is prescribed.
Simulations are performed by considering = (0, 1) over the time domain I = [0, 100), under the following assumptions: • Average constant acceleration (middle point rule) for Newmark scheme, i.e., • Material parameters: κ = ρ E ∈ K = {25, 100, 400, 500}; . In particular, the time taken by the wave to travel along the whole beam (i.e., to propagate from the left-to the rightendpoint) is equal to t = 1/c ∈ {5, 10, 20, 10 √ 5}. Now, since we have set the number of the macro subintervals to M = 5, it turns out that the microscale time variable τ varies in the interval [0, 20). As a consequence, we register two complete wave paths (i.e., two round trips where three reflections take place) for κ = 25, a single complete wave path (with a unique reflection) when κ = 100, and a single one-way path for κ = 400. On the contrary, when κ = 500, the wave propagates too slow with respect to the microscale  adopted, causing a delay between macro-dofs and reflection. This discrepancy is evident, for instance, by the low accuracy characterizing the contour plot in the bottom-right panel of Fig. 7.
These comments find a further confirmation in the trend characterizing the PGD modes associated with the temporal micro-and macro-scales, and with the spatial one, collected in Figs. 8, 9 and 10, respectively for the first three values of κ. In particular, for κ = 25, 100, the information about the periodicity of the solution are already rich enough at the microscale, which allows us to capture the complete wave path (single and double, respectively), whereas the modes over the macroscale are almost constant (only information about wave amplitude are gathered). On the contrary, when κ = 400, the microscale follows only a one-way migration and the information about the wave reflection is recovered at the macro level through more informative macroscale modes.
To sum up, we can state that the choice of the macroscale (and, consequently, of the microscale) is crucial to properly manage the periodicity of the solution. A rough selection of M and of m + 1 may lead to an inaccurate solution, unless a sufficiently large number of temporal modes is adopted for the PGD reconstruction. Convergence of the PGD procedure can be also deteriorated by a not proper choice of parameters M and m + 1. This issue is highlighted by Fig. 11 which shows the trend of the L 2 ( )-norm of the relative error associated with the PGD approximation as a function of the number of the PGD enrichment iterations. As a reference, we adopt the solution to (24)-(25) coinciding with a centered finite difference discretization (in space) combined with the Newmark scheme (in time). In particular, the slow decay of the PGD error in the right panel, for κ = 500, can be ascribed to the decomposition adopted at the macro-and at the micro-scale, which does not match the periodicity of the wave propagation.

Conclusions
The present work successfully accomplished addressing problems exhibiting different time scales within a multi-time separated Proper Generalized Decomposition. The classical one-dimensional time is converted into a multidimensional temporal variable, where different time coordinates correspond to different levels of resolution. This strategy leads to a considerable reduction in terms of computational complexity. Moreover, the issue related to the continuity enforcement between different scales, that motivated several works in the past, is here addressed in a very simple and efficient way, by writing the tensor decomposition in order to reproduce the tensor form of the single time scale discretization. The proposed method is numerically verified first on the one-dimensional heat problem and successively generalized to the wave propagation in elastodynamics.
As a first extension to this work, we are planning to apply the method here developed to visco-elastodynamic equations in a multiscale framework, and, successively, to the nonlinear setting of the visco-elastic-plastic model. In this context, our main interest is for the response of visco-elastic-plastic materials subject to long-term cyclic fatigue loadings [26][27][28] and, in particular, under Very High Cycle Fatigue regimes [29]. Separated space-time representations and hyper-reduction of mechanical models have already been successfully applied in former studies to elastic-plastic problems [30,31], as well as to nonlinear visco-elastic systems [5,[32][33][34]. The idea is to extend the representations provided in this literature to a time multiscale framework, by further tensorization of the time coordinate, following the approach proposed in this paper. Moreover, some theoretical issues characterizing the multi-time PGD method here introduced, such as the impact of nonperiodicity in time caused by irreversible plastic deformations, deserve to be addressed with a more thorough analysis.