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

Angelo Pasquale, Amine Ammar, Antonio Falcó, Simona Perotto, Elías Cueto, Jean-Louis Duval, Francisco Chinesta

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Angelo Pasquale, Amine Ammar, Antonio Falcó, Simona Perotto, Elías Cueto, et al. A separated representation involving multiple time scales within the Proper Generalized Decomposition framework. Advanced Modeling and Simulation in Engineering Sciences, 2021, 8 (26), pp.1-22. 10.1186/s40323-021-00211-7. hal-04096707

HAL Id: hal-04096707

https://hal.science/hal-04096707

Submitted on 29 Jun 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

Open Access

<figure><img src="image_4.png" /></figure>

Check for
updates

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

Angelo Pasquale¹,²*, Amine Ammar¹, Antonio Falcó³, Simona Perotto⁴, Elías Cueto⁵, Jean-Louis Duval⁶ and Francisco Chinesta⁷

**Abstract**

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.

**Keywords:** Model Order Reduction, Time multiscale, Proper Generalized Decomposition, Thermal problems, Elastodynamics

**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.

<figure><img src="image_5.png" /></figure>

<footnote>*Correspondence:
angelo.pasquale@ensam.eu
¹ESI Chair, LAMPA, Arts et Métiers
ParisTech, 2 Boulevard du
Ronceray BP 93525, 49035
Angers Cedex 01, France
Full list of author information is
available at the end of the article</footnote>
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–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,<fnref n="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.

## Time multiscale with PGD

### 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–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 $\Omega$ 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:

$$ \begin{array}{l} \text{Find } \psi \in V(\Omega) \text{ s.t.} \\ a(\psi(x_1, \dots, x_d), \psi^*(x_1, \dots, x_d)) = b(\psi^*(x_1, \dots, x_d)) \quad \forall \psi^* \in V(\Omega), \end{array} \tag{1} $$

with $V(\Omega)$ 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:

$$ {}^{t}\psi^* A \psi = {}^{t}\psi^* B, \tag{2} $$

where the operator on the left-hand side and the right-hand side member are expressed in a separated form as

$$ A = \sum_{k=1}^{N_A} \mathbf{A}_1^k \otimes \mathbf{A}_2^k \otimes \dots \otimes \mathbf{A}_d^k, \quad B = \sum_{k=1}^{N_B} \mathbf{b}_1^k \otimes \mathbf{b}_2^k \otimes \dots \otimes \mathbf{b}_d^k. \tag{3} $$

The approximated PGD solution $\psi$ is sought in the discrete separated form

$$ \psi = \sum_{j=1}^{N_F} \alpha_j \mathbf{f}_1^j \otimes \mathbf{f}_2^j \otimes \dots \otimes \mathbf{f}_d^j \tag{4} $$

and an analogous representation is adopted for the test function $\psi^*$, being

$$ \psi^* = \mathbf{f}_1^i \otimes \mathbf{f}_2^i \otimes \dots \otimes \mathbf{f}_d^i, $$

with $i = 1, \dots, N_F$.

Solution $\psi$ 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

<footnote n="1">¹Even though we only focus on tensor-based methods, we refer to the recent textbooks [10–12] for the latest developments in Model Order Reduction with several applications to multiscale and multiphysics problems.</footnote>
$\epsilon = ||\mathcal{A} \psi - \mathcal{B}||_2$ is small enough, where $||\cdot||_2$ denotes the standard Frobenius norm of a tensor of order $d$.

The projection stage consists of finding the set of coefficients $\alpha_j$ in (4) verifying the relations

$$ \sum_{j=1}^{N_F} H_{ij} \alpha_j = J_i \quad \forall i = 1, \dots, N_F \qquad (5) $$

where

$$ H_{ij} = \sum_{k=1}^{N_A} (\mathbf{t}\mathbf{f}_1^i \mathbf{A}_1^k \mathbf{f}_1^j) \cdot (\mathbf{t}\mathbf{f}_2^i \mathbf{A}_2^k \mathbf{f}_2^j) \dots (\mathbf{t}\mathbf{f}_d^i \mathbf{A}_d^k \mathbf{f}_d^j), $$

$$ J_i = \sum_{k=1}^{N_B} (\mathbf{t}\mathbf{f}_1^i \mathbf{b}_1^k) \cdot (\mathbf{t}\mathbf{f}_2^i \mathbf{b}_2^k) \dots (\mathbf{t}\mathbf{f}_d^i \mathbf{b}_d^k). $$

The enrichment stage includes new candidates for enriching the reduced separated approximation basis, so that $\psi$ can be updated as

$$ \psi = \sum_{j=1}^{N_F} \alpha_j \mathbf{f}_1^j \otimes \mathbf{f}_2^j \otimes \dots \otimes \mathbf{f}_d^j + \underbrace{\mathbf{r}_1 \otimes \mathbf{r}_2 \otimes \dots \otimes \mathbf{r}_d}_{=:\psi_{\mathbf{r}}} \qquad (6) $$

and (2) is replaced by the system

$$ \mathbf{t}^{\psi} \mathcal{A} \psi_{\mathbf{r}} + \mathbf{t}^{\psi} \mathcal{A} \psi_{\mathbf{f}} = \mathbf{t}^{\psi} \mathcal{B}. \qquad (7) $$

In particular, within a fixed point alternating direction algorithm, at each iteration we look for the computation of a single discrete function, $\mathbf{r}_j$, all the other components of $\psi_{\mathbf{r}}$ being assumed known, and after setting $\psi^*$ to $\mathbf{r}_1 \otimes \dots \otimes \mathbf{r}_{j-1} \otimes \mathbf{r}_j^* \otimes \mathbf{r}_{j+1} \dots \otimes \mathbf{r}_d$. This strategy leads to solve the linear system

$$ \mathbf{K} \mathbf{r}_j + \mathbf{v} = \tilde{\mathbf{v}}, \qquad (8) $$

where

$$ \mathbf{K} = \sum_{k=1}^{N_A} \left( \mathbf{A}_j^k \prod_{\substack{h=1 \\ h \neq j}}^{d} \mathbf{r}_h \mathbf{A}_h^k \mathbf{r}_h \right), $$

$$ \mathbf{v} = \sum_{i=1}^{N_F} \sum_{k=1}^{N_A} \left( \alpha_i \mathbf{A}_j^k \mathbf{f}_j^i \prod_{\substack{h=1 \\ h \neq j}}^{d} \mathbf{r}_h \mathbf{A}_h^k \mathbf{f}_h^i \right), \quad \tilde{\mathbf{v}} = \sum_{k=1}^{N_B} \left( \mathbf{b}_j^k \prod_{\substack{h=1 \\ h \neq j}}^{d} \mathbf{r}_h \mathbf{b}_h^k \right). $$

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 $\psi$. 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 $\alpha_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:

$$ ^t\psi^{*t}AA\psi = ^t\psi^{*t}AB, $$

and the associated projection stage becomes

$$ \sum_{j=1}^{N_F} H_{ij}^S \alpha_j = J_i^S \quad \forall i = 1, \dots, N_F $$

where

$$ H_{ij}^S = \sum_{k'=1}^{N_A} \sum_{k=1}^{N_A} ({}^t\mathbf{f}_1^{i't} \mathbf{A}_1^{k'} \mathbf{A}_1^k \mathbf{f}_1^j) \cdot ({}^t\mathbf{f}_2^{i't} \mathbf{A}_2^{k'} \mathbf{A}_2^k \mathbf{f}_2^j) \dots ({}^t\mathbf{f}_d^{i't} \mathbf{A}_d^{k'} \mathbf{A}_d^k \mathbf{f}_d^j), $$

$$ J_i^S = \sum_{k'=1}^{N_A} \sum_{k=1}^{N_B} ({}^t\mathbf{f}_1^{i't} \mathbf{A}_1^{k'} \mathbf{b}_1^k) \cdot ({}^t\mathbf{f}_2^{i't} \mathbf{A}_2^{k'} \mathbf{b}_2^k) \dots ({}^t\mathbf{f}_d^{i't} \mathbf{A}_d^{k'} \mathbf{b}_d^k). $$

Analogously, system (7) will be replaced by the new one

$$ ^t\psi^{*t}AA\psi_r + ^t\psi^{*t}AA\psi_f = ^t\psi^{*t}AB, $$

so that system (8) takes the form

$$ \mathbf{K}^S \mathbf{r}_j + \mathbf{v}^S = \tilde{\mathbf{v}}^S, $$

with

$$ \mathbf{K}^S = \sum_{k'=1}^{N_A} \sum_{k=1}^{N_A} \left( {}^t\mathbf{A}_j^{k'} \mathbf{A}_j^k \prod_{\substack{h=1 \\ h \neq j}}^{d} {}^t\mathbf{r}_h {}^t\mathbf{A}_h^{k'} \mathbf{A}_h^k \mathbf{r}_h \right), $$

$$ \mathbf{v}^S = \sum_{i=1}^{N_F} \sum_{k'=1}^{N_A} \sum_{k=1}^{N_A} \left( \alpha_i {}^t\mathbf{A}_j^{k'} \mathbf{A}_j^k \mathbf{f}_j^i \prod_{\substack{h=1 \\ h \neq j}}^{d} {}^t\mathbf{r}_h {}^t\mathbf{A}_h^{k'} \mathbf{A}_h^k \mathbf{f}_h^i \right), $$

$$ \tilde{\mathbf{v}}^S = \sum_{k'=1}^{N_A} \sum_{k=1}^{N_B} \left( {}^t\mathbf{A}_j^{k'} \mathbf{b}_j^k \prod_{\substack{h=1 \\ h \neq j}}^{d} {}^t\mathbf{r}_h {}^t\mathbf{A}_h^{k'} \mathbf{b}_h^k \right). $$

## 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

$$u(x) = \sum_{i=1}^{N} N_i(x) u_i \sum_{j=1}^{J} G_j(\xi(x - x_i)) g_j, \qquad (9)$$

where the collection, $\{x_i\}_{i=1}^N$, 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$, $\xi(x - x_i)$ is a dependent variable characterized by an offset based on $x_i$, and $G_j(\xi)$ 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 $\xi$. 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 $\Omega_x = (0, L)$ over the time window $\Omega_t = (0, T_f)$, and completed with non homogeneous Dirichlet data,

$$\left\{
\begin{array}{ll}
\displaystyle \frac{\partial u}{\partial t} - k \frac{\partial^2 u}{\partial x^2} = f & \text{in } \Omega_x \times \Omega_t \\
\\
u(x, 0) = u_0(x) & \text{in } \Omega_x \\
\\
u(0, t) = u(L, t) = u_g(t) & \text{in } \Omega_t,
\end{array}
\right.
\qquad (10)$$

with $k$ the diffusivity coefficient and $f = f(x, t)$ the source term.

A standard time-marching approach requires a discretization of $\Omega_t$ by means of $N_t$ times. On the contrary, the technique here presented is based on splitting the whole time interval $(0, T_f)$ into $M$ macro subintervals $\Omega_{i+1} = [T_i, T_{i+1}]$, for $i = 0, \dots, M-1$, of constant amplitude $\Delta T$, by choosing a sequence $\{T_i\}_{i=0}^M$ 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, $\tau(\Omega_1)$, is introduced so that the subinterval $\Omega_1$ is discretized by involving $m+1$ micro-dofs, $\{\tau_j^1\}_{j=0}^m$ and setting $\tau_0^1 = T_0$. Also the discretization at the microscale is uniform, with a constant time step set to $\Delta\tau$ throughout $\Omega_1$. Moreover, the last micro-dof is internal to $\Omega_1$, being $\tau_m^1 = T_1 - \Delta\tau$. By means of a tensor product, the micro discretization in $\Omega_1$ is replicated throughout the other subintervals $\Omega_{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 $\Omega_{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).
<figure><img src="image_6.png" /><figcaption>Fig. 1 Construction of macro and micro scales</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 2 Tensor product of the macro and micro dofs</figcaption></figure>

If we discretize the time derivative in (10) by means of a standard finite difference formula, we have:

$$ \forall n \in \{0, 1, \dots, N_t - 1\}, \quad \frac{\partial u}{\partial t} \approx \frac{1}{\Delta t} (u_{n+1} - u_n), \qquad (11) $$

that we can write also as

$$ \frac{\partial u}{\partial t} \approx \mathbf{E}_{N_t} \mathbf{u}, \qquad (12) $$
with

$$ E_{N_t} = \frac{1}{\Delta t} \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & 0 \\ -1 & 1 & 0 & \cdots & 0 & 0 \\ 0 & -1 & 1 & 0 & 0 & \cdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & \cdots & 0 & -1 & 1 & 0 \\ 0 & 0 & \cdots & 0 & -1 & 1 \\ 0 & 0 & \cdots & 0 & 0 & -1 & 1 \end{pmatrix}_{N_t} \quad \text{and} \quad \mathbf{u} = \begin{pmatrix} u_1 \\ u_2 \\ \vdots \\ u_{N_t} \end{pmatrix}, \qquad (13) $$

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

$$ \frac{\partial}{\partial \tau} u(T, \tau), $$

where $T \in \{T_i\}_{i=0}^{M-1}$ is adopted to identify the subinterval $\Omega_{i+1}$ of interest, while $\tau \in [T_0, T_1)$ represents the actual temporal independent variable, spanning $\Omega_1$, i.e, via tensor product, the generic subinterval $\Omega_{i+1}$. This leads to perform the time derivation with respect to variable $\tau$.

With a view to the application of PGD for approximating problem (10) with the two-scale 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$, $\frac{\partial^0}{\partial T^0}$, with the identity operator, so that, $\frac{\partial^0}{\partial T^0}f(T) = f(T)$, and, as a consequence, identity

$$ \frac{\partial}{\partial \tau} u(T, \tau) = \left( \frac{\partial^0}{\partial T^0} \otimes \frac{\partial}{\partial \tau} \right) u(T, \tau) $$

holds true. This equality suggested the authors in [1] to propose the approximation

$$ \frac{\partial}{\partial t} \approx \frac{\partial^0}{\partial T^0} \otimes \frac{\partial}{\partial \tau}, \qquad (14) $$

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

$$ E_{N_t} \approx E_{T^0, \tau} = E_{m+1} \otimes I_M = \begin{pmatrix} E_{m+1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \cdots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & E_{m+1} \end{pmatrix}_M, \qquad (15) $$

where now the subscript $M$ refers to the number of blocks constituting matrix $E_{T^0, \tau}$, while matrix $E_{m+1}$, sized according to discretization at the micro-scale, is defined as in (13), after replacing $\Delta t$ with $\Delta \tau$, and replicated throughout the subintervals $\Omega_{i+1}$ by means of the Kronecker product with the discrete identity operator $I_M$.<fnref n="2" />

$$ \mathbf{A} \otimes \mathbf{B} = \begin{pmatrix} B_{1,1} & \mathbf{A} & B_{1,2} & \mathbf{A} & \cdots & B_{1,n} & \mathbf{A} \\ B_{2,1} & \mathbf{A} & B_{2,2} & \mathbf{A} & \cdots & B_{2,n} & \mathbf{A} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ B_{m,1} & \mathbf{A} & B_{m,2} & \mathbf{A} & \cdots & B_{m,n} & \mathbf{A} \end{pmatrix} \quad \text{and} \quad \mathbf{a} \otimes \mathbf{b} = \begin{pmatrix} b_1 \mathbf{a} \\ \vdots \\ b_m \mathbf{a} \end{pmatrix}, $$

<footnote n="2">²We observe that symbol $\otimes$ refers to a revision of the standard Kronecker product, being meant</footnote>
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, $\Omega_{i+1}$, to the first dof associated with the macro-scale of the subinterval $\Omega_{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

$$ E_{N_t} = \frac{1}{\Delta\tau} \left( \begin{array}{c} \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ -1 & 1 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & -1 & 1 \\ 0 & \cdots & 0 & -1 & 1 \end{pmatrix}_{m+1} \\ + \frac{1}{\Delta\tau} \begin{pmatrix} 0 & 0 & 0 & \cdots & 0 & -1 \\ 0 & 0 & 0 & 0 & \cdots & 0 \\ 0 & \ddots & \ddots & \ddots & 0 & \vdots \\ \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 & 0 \end{pmatrix}_{m+1} \end{array} \right) \otimes \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ \vdots & \ddots & 1 & 0 \\ 0 & \cdots & 0 & 1 \\ 0 & \cdots & 0 & 1 & 0 \end{pmatrix}_M \tag{16} $$

which, in a compact form, reads as

$$ E_{N_t} = E_{m+1} \otimes I_M + C_{m+1} \otimes L_M \tag{17} $$

where matrices $C_{m+1}$ and $L_M$ are defined as

$$ (C_{m+1})_{ij} = -\frac{1}{\Delta\tau} \begin{cases} 1 & \text{if } i=1, j=m+1 \\ 0 & \text{otherwise,} \end{cases} \quad (L_M)_{ij} = \delta_{ij+1}, \tag{18} $$

respectively with $\delta_{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>Footnote 3 continued</footnote>

<footnote>respectively. This is a common choice in the PGD community to preserve the same order of the variables in the separated form [23].</footnote>
Coming back to the approximation of the heat equation in (10), a classical discrete PGD approximation, based on a space-time separated form, reads:

$$ \text{Find } \mathbf{U}_H = \sum_{k=1}^{N_U} \mathbf{u}_1^k \otimes \mathbf{u}_2^k \text{ such that } \mathcal{A}_H \mathbf{U}_H = \mathcal{B}_H, \qquad (19) $$

where $N_U$ denotes the number of the PGD modes, index 1 refers to the space dependence, index 2 is associated with the time, $\mathcal{B}_H$ is a suitable discrete separated approximation of the source term $f(t, x)$ evaluated at some discrete points, $\mathcal{A}_H$ is a discrete separated representation of the differential operators in (10), built by exploiting the tensorial identity

$$ \frac{\partial}{\partial t} - k \frac{\partial^2}{\partial x^2} = \left( \frac{\partial}{\partial t} \otimes \frac{\partial^0}{\partial x^0} \right) - k \left( \frac{\partial^0}{\partial t^0} \otimes \frac{\partial^2}{\partial x^2} \right), \qquad (20) $$

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 $\mathcal{A}_H$ can be separated as

$$ \mathcal{A}_H = \sum_{k=1}^{2} \mathbf{A}_1^k \otimes \mathbf{A}_2^k \text{ with } \mathbf{A}_1^1 = \mathbf{I}_{N_x}, \mathbf{A}_2^1 = \mathbf{E}_{N_t}, \mathbf{A}_1^2 = -k \mathbf{D}_{N_x}, \mathbf{A}_2^2 = \mathbf{I}_{N_t}, \qquad (21) $$

where $\mathbf{E}_{N_t}$ denotes the matrix associated with the implicit Euler scheme, $\mathbf{I}_{N_x}$ and $\mathbf{I}_{N_t}$ are the identity matrices in space and time, respectively and $\mathbf{D}_{N_x}$ is the matrix of the discrete Laplacian. Figure 3 shows the sparsity pattern associated with four matrices $\mathbf{A}_i^k$.

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

$$ \frac{\partial}{\partial t} - k \frac{\partial^2}{\partial x^2} = \left( \frac{\partial}{\partial t} \otimes \frac{\partial^0}{\partial x^0} \right) - k \left( \frac{\partial^0}{\partial t^0} \otimes \frac{\partial^2}{\partial x^2} \right) \\ \approx \left( \frac{\partial^0}{\partial T^0} \otimes \frac{\partial}{\partial \tau} \otimes \frac{\partial^0}{\partial x^0} \right) - k \left( \frac{\partial^0}{\partial T^0} \otimes \frac{\partial^0}{\partial \tau^0} \otimes \frac{\partial^2}{\partial x^2} \right), $$

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 $\mathcal{A}_H$ in (19) can be separated as

$$ \mathcal{A}_H = \sum_{k=1}^{3} \mathbf{A}_1^k \otimes \mathbf{A}_2^k \otimes \mathbf{A}_3^k, \qquad (22) $$

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

$$ \mathbf{A}_1^1 = \mathbf{A}_1^2 = \mathbf{I}_{N_x}, \quad \mathbf{A}_1^3 = -k \mathbf{D}_{N_x}, $$
<figure><img src="image_8.png" /><figcaption>Fig. 3 Sparsity pattern of discrete PGD operators associated with the separated representation (21), for $N_x = 101$ and $N_t = 1000$</figcaption></figure>

whereas the two-scale time discretization leads us to adopt the modified tensorial decomposition in (16),

$$E_{N_t} = A_2^1 \otimes A_3^1 + A_2^2 \otimes A_3^2,$$

with $A_2^1 = E_{m+1}$ the matrix characterizing the implicit Euler discretization at the microscale; $A_3^1 = A_3^3 = I_M$ the discrete identity matrix associated with the subintervals $\Omega_{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_3^2 = L_M$ the lower shift matrix, sized according to the macroscale, with the first subdiagonal of elements all equal to one; $A_2^3 = 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
<figure><img src="image_9.png" /><figcaption>**Fig. 4** Sparsity pattern of the discrete PGD operators associated with the separated representation (22), for $N_x = 101$, $m + 1 = 100$ and $M = 10$</figcaption></figure>

about equal to $M^{-1}$, for $M \ll m + 1$. Notice that the time decomposition in Fig. 4 suits the time grid in Fig. 3, being $M(m + 1) = N_t = 1000$.

## 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 one-dimensional spatial domain which models, for instance, the traction-compression waves travelling along a linear elastic medium, $\Omega = (0, L)$, during time interval $I = [0, T)$. The scalar displacement field is denoted by $u = u(x, t)$, for $x \in \Omega$ and $t \in 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, $\partial\Omega_u \times I$ and $\partial\Omega_f \times I$, of the boundary of the space-time cylinder, respectively, with $\partial\Omega_f \cup \partial\Omega_u = \partial\Omega$ and $\partial\dot{\Omega}_f \cap \partial\dot{\Omega}_u = \emptyset$ (notice that, in the specific context at hand, $\partial\Omega_u$ and $\partial\Omega_f$ reduce to the endpoints of the domain $\Omega$, 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, $\rho$, an elasticity modulus, $E$, and a section, $A$. Thus, the strong formulation of the considered transient elastodynamic problem reads: find $u : \Omega \times I \to \mathbb{R}$, such that

$$ 
\left\{
\begin{array}{ll}
\rho \frac{\partial^2 u}{\partial t^2} = E \frac{\partial^2 u}{\partial x^2} & \text{in } \Omega \times I, \\
u = u^d & \text{on } \partial \Omega_u \times I, \\
\frac{\partial u}{\partial x} = \frac{f^d}{EA} & \text{on } \partial \Omega_f \times I, \\
u(x, 0) = u_0(x) & \text{on } \Omega \times \{0\}, \\
\frac{\partial u}{\partial t}(x, 0)(x, 0) = v_0(x) & \text{on } \Omega \times \{0\}.
\end{array}
\right.
\qquad (23) $$

To simplify the notation, in the following we adopt the convention $\dot{u}$ and $\ddot{u}$ to denote the first- and the second-order time derivative, $\frac{\partial u}{\partial t}$ and $\frac{\partial^2 u}{\partial t^2}$, of $u$ respectively.

We employ a standard centered finite difference scheme to approximate the space, after discretizing $\Omega$ with $N_x$ dofs, so that the associated semi-discrete formulation is: find $\mathbf{u} : I \to \mathbb{R}^{N_x}$, such that

$$ 
\left\{
\begin{array}{l}
\mathbf{M} \ddot{\mathbf{u}}(t) + \mathbf{K} \mathbf{u}(t) = \mathbf{F}(t), \\
\mathbf{u}(0) = \mathbf{u}_0, \\
\dot{\mathbf{u}}(0) = \mathbf{v}_0,
\end{array}
\right.
\qquad (24) $$

where $\mathbf{M} = \frac{\rho}{E} \mathbf{I}_{N_x}$ and $\mathbf{K} = - \mathbf{D}_{N_x}$ denote the mass and the stiffness matrix, respectively, with $\mathbf{I}_{N_x}$ and $\mathbf{D}_{N_x}$ the identity and discrete Laplacian operators, $\mathbf{u}_0$ and $\mathbf{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 $\mathbf{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} = \Delta t$, for $i = 1, \dots, N_t$. Following [23], we apply the Newmark scheme to (24), so that we have to solve the fully discretized problem: find $\mathbf{u} \in \mathbb{R}^{N_x \times N_t}$, such that

$$ (\mathbf{K} \otimes \mathbf{N}_1 + \mathbf{M} \otimes \mathbf{N}_2) : \mathbf{u} = (\mathbf{I}_{N_x} \otimes \mathbf{N}_1) : \mathbf{F} + (\mathbf{M} \mathbf{u}_0) \otimes \mathbf{N}_3 + (\mathbf{M} \mathbf{v}_0) \otimes \mathbf{N}_4 \\ + (\mathbf{F}(t_0) - \mathbf{K} \mathbf{u}_0) \otimes \mathbf{N}_5, \qquad (25) $$

where : denotes the two times contracted product between tensors<fnref n="3" />,

$$ 
\begin{gathered}
\mathbf{N}_1 = \begin{pmatrix}
\beta & 0 & \cdots & 0 \\
a & \ddots & \ddots & \vdots \\
b & \ddots & \ddots & \vdots \\
0 & \ddots & \ddots & 0 \\
\cdots & b & a & \beta
\end{pmatrix}_{N_t}, \qquad \mathbf{N}_2 = \frac{1}{\Delta t^2} \begin{pmatrix}
1 & 0 & \cdots & 0 \\
-2 & \ddots & \ddots & \vdots \\
1 & \ddots & \ddots & \vdots \\
0 & \ddots & \ddots & 0 \\
\cdots & 0 & 1 & -2
\end{pmatrix}_{N_t}, \\
\mathbf{N}_3 = \frac{1}{\Delta t^2} \begin{pmatrix}
1 \\
-1 \\
0 \\
\vdots \\
0
\end{pmatrix}_{N_t},
\end{gathered}
$$

<footnote n="3">³The two times contracted product : between second or fourth order tensors is defined as $\mathbf{A} : \mathbf{B} = \sum_{i,j} A_{ij} B_{ij}$.</footnote>
$$ N_4 = \frac{1}{\Delta t} \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}_{N_t}, \quad N_5 = \begin{pmatrix} \frac{1}{2} - \beta \\ \frac{1}{2} + \beta - \gamma \\ 0 \\ \vdots \\ 0 \end{pmatrix}_{N_t} \qquad (26) $$

with $0 \le \gamma \le 1$ and $0 \le \beta \le \frac{1}{2}$ the parameters involved by the Newmark algorithm, $a = 1/2 - 2\beta + \gamma$, $b = 1/2 + \beta - \gamma$, $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, \dots, 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

$$ \text{Find } \mathbf{U}_E = \sum_{k=1}^{N_{II}} \mathbf{u}_1^k \otimes \mathbf{u}_2^k \text{ such that } \mathcal{A}_E \mathbf{U}_E = \mathcal{B}_E, \qquad (27) $$

where

$$ \mathcal{A}_E = \sum_{k=1}^{2} \mathbf{A}_1^k \otimes \mathbf{A}_2^k, \quad \text{with} \quad \mathbf{A}_1^1 = \mathbf{M}, \quad \mathbf{A}_2^1 = \mathbf{N}_2, \quad \mathbf{A}_1^2 = \mathbf{K}, \quad \mathbf{A}_2^2 = \mathbf{N}_1, \qquad (28) $$

the same notation as in (19) being adopted for the subscripts identifying the space and the time, respectively, while $\mathcal{B}_E$ simplifies to

$$ \mathcal{B}_E = \mathbf{b}_1 \otimes \mathbf{b}_2, \quad \text{with} \quad \mathbf{b}_1 = \frac{\rho}{E} \mathbf{v}_0, \quad \mathbf{b}_2 = \mathbf{N}_4 $$

after assuming $\mathbf{F} = \mathbf{0}$ (i.e., no external forces) and $\mathbf{u}_0 = \mathbf{0}$ (i.e., a vanishing initial displacement).

The discrete operators involved in the definition of $\mathcal{A}_E$ have the sparsity pattern shown in Fig. 5, when choosing $\Omega = (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, $\mathbf{E}_{N_t}$, is here extended to the Newmark discrete operators, $\mathbf{N}_1$ and $\mathbf{N}_2$, in (26), so that it holds

$$ \mathbf{N}_1 = \begin{pmatrix} \beta & 0 & 0 & \dots & 0 \\ a & \beta & 0 & 0 & 0 \\ b & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & b & a & \beta \\ 0 & \dots & b & a & \beta \end{pmatrix}_{m+1} \otimes \begin{pmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \vdots \\ \vdots & \ddots & 1 & 0 \\ 0 & \dots & 0 & 1 \end{pmatrix}_M $$
$$ + \begin{pmatrix} 0 & 0 & \cdots & 0 & b & a \\ 0 & 0 & \ddots & 0 & \ddots & b \\ 0 & \ddots & \ddots & \ddots & 0 & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 & 0 \end{pmatrix}_{m+1} \otimes \begin{pmatrix} 0 & 0 & \cdots & 0 \\ 1 & 0 & \ddots & \vdots \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 1 & 0 \\ 0 & \cdots & 0 & 1 & 0 \\ 0 & \cdots & 0 & 0 & 1 & 0 \end{pmatrix}_{M}, \quad (29) $$

$$ \mathbf{N}_2 = \frac{1}{\Delta \tau^2} \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 \\ -2 & 1 & 0 & \cdots & 0 \\ 1 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 1 & -2 & 1 & 0 \\ 0 & \cdots & 1 & -2 & 1 & 0 \end{pmatrix}_{m+1} \otimes \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots \\ \vdots & \ddots & 1 & 0 \\ 0 & \cdots & 0 & 1 & 0 \\ 0 & \cdots & 0 & 0 & 0 & 0 \end{pmatrix}_{M} \\ + \frac{1}{\Delta \tau^2} \begin{pmatrix} 0 & 0 & \cdots & 0 & 1 & -2 \\ 0 & 0 & \ddots & 0 & \ddots & 1 \\ 0 & \ddots & \ddots & \ddots & 0 & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & 0 & 0 & 0 \\ 0 & 0 & \cdots & 0 & 0 & 0 \end{pmatrix}_{m+1} \otimes \begin{pmatrix} 0 & 0 & \cdots & 0 \\ 1 & 0 & \ddots & \vdots \\ 0 & 1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 1 & 0 \\ 0 & \cdots & 0 & 1 & 0 \\ 0 & \cdots & 0 & 0 & 1 & 0 \end{pmatrix}_{M}, \quad (30) $$

$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

$$ \mathbf{N}_l = \mathbf{N}_{l,m+1} \otimes \mathbf{I}_M + \mathbf{C}_{l,m+1} \otimes \mathbf{L}_M \quad \text{with } l = 1, 2, \quad (31) $$

where matrices $\mathbf{N}_{1,m+1}$ and $\mathbf{N}_{2,m+1}$ are defined exactly as $\mathbf{N}_1$ and $\mathbf{N}_2$, respectively except for the size (equal to $m+1$ instead of $N_t$); $\mathbf{I}_M$ is the discrete identity matrix associated with the macroscale; $\mathbf{L}_M$ is the lower shift matrix defined as in (18); the correction matrices $\mathbf{C}_{1,m+1}$ and $\mathbf{C}_{2,m+1}$ are given by

$$ (\mathbf{C}_{1,m+1})_{ij} = \begin{cases} a & \text{if } i=1, j=m+1 \\ b & \text{if } i=1, j=m \text{ or } i=2, j=m+1, \\ 0 & \text{otherwise,} \end{cases} $$

$$ (\mathbf{C}_{2,m+1})_{ij} = \frac{1}{\Delta \tau^2} \begin{cases} -2 & \text{if } i=1, j=m+1 \\ 1 & \text{if } i=1, j=m \text{ or } i=2, j=m+1, \\ 0 & \text{otherwise.} \end{cases} $$

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 $\mathcal{A}_E$ in (27), after separating this operator in terms of the coordinates $x$, $T$ and $\tau$, and by exploiting decompositions (31), so that we have

$$ \mathcal{A}_E = \sum_{k=1}^{4} \mathbf{A}_1^k \otimes \mathbf{A}_2^k \otimes \mathbf{A}_3^k, \quad (32) $$
<figure><img src="image_10.png" /><figcaption>Fig. 5 Sparsity pattern of discrete PGD operators associated with the separated representation (28), for $N_x = 51$ and $N_t = 2000$</figcaption></figure>

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

$$A_1^1 = A_1^2 = M = \frac{\rho}{E} I_{N_x}, \quad A_1^3 = A_1^4 = K = -D_{N_x}.$$

The operators associated with the time microscale are

$$A_2^1 = N_{2,m+1}, \quad A_2^2 = C_{2,m+1}, \quad A_2^3 = N_{1,m+1}, \quad A_2^4 = C_{1,m+1},$$

while the ones related to the temporal macroscale are

$$A_3^1 = A_3^3 = I_M, \quad A_3^2 = A_3^4 = L_M.$$

In particular, the time microscale is characterized by the Toeplitz matrices, $A_2^1$ and $A_2^3$, and by the correction matrices, $A_2^2$ and $A_2^4$, with only three not null elements; the diagonal matrices $A_3^1$ and $A_3^3$ and the lower shift matrices $A_2^2$ and $A_2^4$ identify the time macroscale. The sparsity pattern associated with the 12 operators in (32) is shown in Fig. 6, for $\Omega = (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

$$B_E = \mathbf{b}_1 \otimes \mathbf{b}_2 \otimes \mathbf{b}_3, \quad \text{with} \quad \mathbf{b}_1 = \frac{\rho}{E} \mathbf{v}_0, \quad \mathbf{b}_2 = \frac{1}{\Delta \tau} \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}_{m+1}, \quad \mathbf{b}_3 = \begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix}_M.$$
<figure><img src="image_11.png" /><figcaption>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$</figcaption></figure>

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_1^2$ in Fig. 5 being preserved by matrices $A_1^1-A_1^2$ and $A_1^3-A_1^4$ 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 Fig. 5, being $N_t = M(m + 1) = 2000$).

## 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 $\Omega = (0, 1)$ over the time domain $I = [0, 100)$, under the following assumptions:

* Average constant acceleration (middle point rule) for Newmark scheme, i.e.,

$$ \beta = \frac{1}{4}, \quad \gamma = \frac{1}{2} \implies a = \frac{1}{2}, \quad b = \frac{1}{4}; $$

* Material parameters: $\kappa = \frac{\rho}{E} \in K = \{25, 100, 400, 500\}$;
<figure><img src="image_12.png" /><figcaption>Fig. 7 Space-time solutions for $κ = 25, 100, 400, 500$ (left-right, top-bottom)</figcaption></figure>

* Sinusoidal initial velocity:

$$v_0(x) = \begin{cases} \frac{1}{100} \sin \left( \frac{\pi x}{0.3} \right)^2 & \text{for } x < 0.3 \\ 0 & \text{otherwise;} \end{cases}$$

* Discretization parameters: $N_x = 51, m + 1 = 400, M = 5$.

Figure 7 collects the space-time contour plot for the propagation of the elastic wave, for each $κ ∈ K$. The values in K have been selected in order to obtain solutions with a different level of detail in terms of behaviour at the microscale. Indeed, quantity $κ$ can be easily related to the speed of the stress wave along the beam, after rewriting the differential equation in (23) as

$$\frac{1}{c^2} \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2},$$

with $c = \sqrt{E/\rho} = \sqrt{1/\kappa} \in C = \{\frac{1}{5}, \frac{1}{10}, \frac{1}{20}, \frac{1}{10\sqrt{5}}\}$. In particular, the time taken by the wave to travel along the whole beam (i.e., to propagate from the left- to the right-endpoint) is equal to $\Delta t = 1/c \in \{5, 10, 20, 10\sqrt{5}\}$. Now, since we have set the number of the macro subintervals to $M = 5$, it turns out that the microscale time variable $\tau$ 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
<figure><img src="image_13.png" /><figcaption>Fig. 8 PGD modes—$κ = 25$</figcaption></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 9 PGD modes—$κ = 100$</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 10 PGD modes—$κ = 400$</figcaption></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 11 Error decay for successive PGD enrichments for $κ = 100$ (left) and $κ = 500$ (right)</figcaption></figure>

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
$\kappa = 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 $\kappa = 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(\Omega)$-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 $\kappa = 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 non-linear 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–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 non-linear visco-elastic systems [5,32–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 non-periodicity in time caused by irreversible plastic deformations, deserve to be addressed with a more thorough analysis.
## Authors' contributions

All the authors participated in the definition of techniques and algorithms. All authors read and approved the final manuscript.

## Funding

The authors acknowledge the support of the ESI Group through its research chair at ENSAM ParisTech and through the project 2019-0060 "Simulated Reality" at the University of Zaragoza. The support of the Spanish Ministry of Science through Grant number CICYT-DPI2017-85139-C2-1-R and by the Regional Government of Aragon and the European Social Fund, are also gratefully acknowledged. This work has been also partially supported by European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, Grant agreement 872442 (ARIA) and by the research project GNCS-INdAM 2020.

## Availability of data and materials

Interested reader can contact authors.

## Declarations

### Competing interests

The authors declare that they have no competing interests.

### Author details

¹ESI Chair, LAMPA, Arts et Métiers ParisTech, 2 Boulevard du Ronceray BP 93525, 49035 Angers Cedex 01, France,

²Politecnico di Milano, Piazza Leonardo da Vinci, 32, 20133 Milano, Italy, ³ESI Chair, Universidad Cardenal Herrera-CEU, San Bartolomé 55, 46115 Alfara del Patriarca (Valencia), Spain, ⁴MOX, Dept. of Mathematics, Politecnico di Milano, Piazza Leonardo da Vinci, 32, 20133 Milano, Italy, ⁵Aragon Institute of Engineering Research, Universidad de Zaragoza, Calle de Mariano Esquillor, s/n, 50018 Zaragoza, Spain, ⁶ESI Group, Parc Icade, Immeuble le Seville, 3 bis, Saarinen, CP 50229, 94528 Rungis Cedex, France, ⁷ESI Group Chair, PIMM, Arts et Métiers ParisTech, 151 Boulevard de l'Hopital, 75013 Paris, France.

Received: 12 April 2021 Accepted: 26 October 2021

Published online: 26 November 2021

## References

1. Ammar A, Chinesta F, Cueto E, Doblare M. Proper generalized decomposition of time-multiscale models. Int J Numer Methods Eng. 2012;90:569-96.

2. Néron D, Ladevèze P. Proper generalized decomposition for multiscale and multiphysics problems. Arch Comput Methods Eng. 2010;17(4):351-72.

3. Passieux JC, Ladevèze P, Néron D. A scalable time-space multiscale domain decomposition method: adaptive time scale separation. Comput Mech. 2010;46(4):621-33.

4. Dureisseix D, Ladevèze P, Schrefler B. A LATIN computational strategy for multiphysics problems-application to poroelasticity. Int J Numer Methods Eng. 2003;56:1489-510.

5. Fritzen F, Hassani MR. Space-time model order reduction for nonlinear viscoelastic systems subjected to long-term loading. Meccanica. 2018;53:1333-55.

6. Maday Y, Turinici G. A parareal in time procedure for the control of partial differential equations. C R Math. 2002;335(4):387-92.

7. Maday Y, Riahi, M, Salomon J. Parareal in time intermediate targets methods for optimal control problems. In: Control and optimization with PDE constraints. Internat Ser Numer Math. Basel: Springer. 2013. p. 79-92.

8. Badias A, Gonzalez D, Alfaro I, Chinesta F, Cueto E. Local proper generalized decomposition. Int J Numer Methods Eng. 2017;112:1715-32.

9. Ibáñez R, Ammar A, Cueto E, Huerta A, Duval JL, Chinesta F. Multiscale proper generalized decomposition based on the partition of unity. Int J Numer Methods Eng. 2019;120(6):727-47.

10. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. Model order reduction applications. Berlin: De Gruyter; 2020.

11. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. Model order reduction: snapshot-based methods and algorithms. Berlin: De Gruyter; 2020.

12. Benner P, Schilders W, Grivet-Talocia S, Quarteroni A, Rozza G, Miguel Silveira L. Model order reduction: system- and data-driven methods and algorithms. Berlin: De Gruyter; 2020.

13. Ammar A. The proper generalized decomposition: a powerful tool for model reduction. Int J Mater Form. 2010;3:89-102.

14. Ammar A, Chinesta F, Díez P, Huerta A. An error estimator for separated representations of highly multidimensional models. Comput Methods Appl Mech Eng. 2010;199:1872-80.

15. Chinesta F, Keunings R, Leygue A. The proper generalized decomposition for advanced numerical simulations. A primer. Springer briefs in applied sciences and technology. Cham: Springer; 2014.

16. Chinesta F, Ammar A, Cueto E. Recent advances and new challenges in the use of the proper generalized decomposition for solving multidimensional models. Arch Comput Methods Eng. 2010;17(4):327-50.

17. Chinesta F, Cueto E. PGD-based modeling of materials. Structures and processes. ESAFORM Bookseries on material forming. Cham: Springer; 2014.

18. Perotto S, Carlino MG, Ballarin F. Model reduction by separation of variables: a comparison between hierarchical model reduction and proper generalized decomposition. In: Sherwin SJ, Moxey D, Peirò J, Vincent PE, Schwab C,
editors. Spectral and high order methods for partial differential equations ICOSAHOM 2018, vol. 134. Series: Lect. Notes Comput. Sci. Eng. Cham: Springer; 2020. p. 61–77.

19. Chinesta F, Ammar A, Lemarchand F, Beauchene P, Boust F. Alleviating mesh constraints: model reduction, parallel time integration and high resolution homogenization. Comput Methods Appl Mech Eng. 2008;197(5):400–13.

20. Chinesta F, Ladevèze P, Cueto E. A short review in model order reduction based on proper generalized decomposition. Arch Comput Methods Eng. 2011;18(4):395–404.

21. Ammar A, Chinesta F, Falcó A. On the convergence of a greedy rank-one update algorithm for a class of linear systems. Arch Comput Methods Eng. 2010;17(4):473–86.

22. de Silva V, Lim LH. Tensor rank and the ill-posedness of the best low-rank approximation problem. SIAM J Matrix Anal Appl. 2006;30(3):1084–127.

23. Boucinha L, Gravouil A, Ammar A. Space-time proper generalized decompositions for the resolution of transient elastodynamic models. Comput Methods Appl Mech Eng. 2013;255:67–88.

24. Lathauwer L, De Moor B. A multi-linear singular value decomposition. J Soc Ind Appl Math. 2000;21:1253–78.

25. Newmark NM. A method of computation for structural dynamics. J Eng Mech. 1959;85:67–94.

26. Oskay C, Fish J. Multiscale modeling of fatigue for ductile materials. Int J Multiscale Comp Eng. 2004;2:329–54.

27. Oskay C, Fish J. Fatigue life prediction using 2-scale temporal asymptotic homogenization. Int J Numer Methods Eng. 2004;61:329–59.

28. Cognard J-Y, Ladevèze P. A large time increment approach for cyclic viscoplasticity. Int J Plast. 1993;9(3):141–57.

29. Furuya Y, Hirukawa H, Takeuchi E. Gigacycle fatigue in high strength steels. Sci Technol Adv Mater. 2019;20(1):643–56.

30. Bergheau J-M, Zuchiatti S, Roux J-C, Feulvarch É, Tissot S, Perrin G. The proper generalized decomposition as a space-time integrator for elastoplastic problems. Meccanica. 2016;344:759–68.

31. Nasri MA, Aguado JV, Ammar A, Cueto E, Chinesta F, Morel F, Robert C, El Arem S. Separated representations of incremental elastoplastic simulations. Key Eng Mater. 2015;651–653:1285–93.

32. Ryckelynck D, Benziane D. Hyper-reduction framework for model calibration in plasticity-induced fatigue. Adv Model Simul Eng Sci. 2016;3:10.

33. Cremonesi M, Néron D, Guidault P, Ladevèze P. A PGD-based homogenization technique for the resolution of nonlinear multiscale problems. Comput Methods Appl Mech Eng. 2013;267:275–92.

34. Ladevèze P, Nouy A. On a multiscale computational strategy with time and space homogenization for structural mechanics. Comput Methods Appl Mech Eng. 2003;192(28–30):3061–87.

**Publisher's Note**

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Submit your manuscript to a SpringerOpen® journal and benefit from:**

► Convenient online submission

► Rigorous peer review

► Open access: articles freely available online

► High visibility within the field

► Retaining the copyright to your article

Submit your next manuscript at ► springeropen.com