# Proper Generalized Decomposition of Multiscale Models

Francisco Chinesta, Amine Ammar, Elías Cueto

► To cite this version:

Francisco Chinesta, Amine Ammar, Elías Cueto. Proper Generalized Decomposition of Multiscale Models. International Journal for Numerical Methods in Engineering, 2010, 83 (8-9), pp.1114-1132. 10.1002/nme.2794 . hal-01007222

HAL Id: hal-01007222

https://hal.science/hal-01007222

Submitted on 30 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

Public Domain
# Proper Generalized Decomposition of Multiscale Models

F. Chinesta¹, A. Ammar², E. Cueto³

¹ EADS Corporate International Chair
GeM – Institut de Recherche en Génie Civil et Mécanique
UMR CNRS - ECN - Université de Nantes
1 rue de la Noe, BP 92101. F-44321 Nantes Cedex 3, France
Francisco.Chinesta@ec-nantes.fr

² Laboratoire de Rhéologie, UMR CNRS - INPG - UJF
1301 rue de la piscine, BP 53 Domaine Universitaire
F-38041 Grenoble Cedex 9, France
Amine.Ammar@ujf-grenoble.fr

³ Group of Structural Mechanics and Material Modelling.
Aragón Institute of Engineering Research, I3A. Universidad de Zaragoza.
Edificio Betancourt. Mara de Luna, 7. E-50012 Zaragoza, Spain
ecueto@unizar.es

**Keywords:** Proper Generalized Decompositions, Separated representations, Finite sums decomposition, Multidimensional models, Model reduction.

**Abstract:** In this paper the coupling of a parabolic model with a system of local kinetic equations is analyzed. A space-time separated representation is proposed for the global model (this is simply the radial approximation proposed by Pierre Ladeveze in the LATIN framework [12]). The main originality of the present work concerns the treatment of the local problem, that is firstly globalized (in space and time) and then fully globalized by introducing a new coordinate related to the different species involved in the kinetic model. Thanks to the non-incremental nature of both discrete descriptions (the local and the global one) the coupling is quite simple and no special difficulties are encountered by using heterogeneous time integrations.
# 1 Introduction

Some models encountered in science and engineering are sometimes defined in multidimensional spaces (e.g. the ones involved in quantum mechanics or kinetic theory descriptions of materials, including complex fluids) that exhibit the redoubtable curse of dimensionality when usual mesh-based discretization techniques are applied. Other times, models involve transient fields, that even when they are defined in three-dimensional physical spaces, they must be solved in large time intervals using very small time steps. Moreover, standard models can become multidimensional if some of the parameters that they involve are considered as new coordinates. This possibility is specially attractive when these coefficients are not well known or they have a stochastic nature.

The difficulty related to the solution of multidimensional models is quite obvious and it needs the proposal of new appropriate strategies able to circumvent the curse of dimensionality. One possibility lies in the use of sparse grids [7]. However, as argued in [1], the use of sparse grid is restricted to models with moderate multidimensionality (up to 20). Another technique able to circumvent, or at least alleviate, the curse of dimensionality consists of using a separated representation of the unknown field (see [20] [6] for some numerical elements on this topic). Basically, the separated representation of a generic function $u(x_1, \dots, x_D)$ (also known as finite sums decomposition) writes:

$$u(x_1, \dots, x_D) \approx \sum_{i=1}^{i=N} F_i^1(x_1) \times \dots \times F_i^D(x_D) \quad (1)$$

This kind of representation is not new, it was widely employed in the last decades in the framework of quantum chemistry. In particular the Hartree-Fock (that involves a single product of functions) and post-Hartree-Fock approaches (as the MCSCF that involves a finite number of sums) made use of a separated representation of the wavefunction [8] [11].

We proposed recently a technique able to construct, in a way completely transparent for the user, the separated representation of the unknown field involved in a multidimensional partial differential equation. This technique, originally described and applied to multi-bead-spring FENE models of polymeric systems in [3], was extended to transient models of such complex fluids in [4]. Other more complex models (involving different couplings and non-linearities) based on the reptation theory of polymeric liquids were analyzed in [15].

Now, we come back to models defined in spaces of moderate dimension ($dD$, $d = 1, 2, 3$) but whose solutions evolve in large time intervals. In this context, if one uses standard incremental time-discretizations, in the general case (models involving time-dependent parameters, non-linear models, ...), one must solve at least a linear system at each time step. When the time step becomes too small as a consequence of the stability requirements, and the simulation time interval is large enough, the simulation becomes inefficient. To illustrate this scenario, one could imagine the simple reaction-diffusion model that describes the degradation of plastic materials, where the characteristic time of the chem-
ical reaction involved in the material degradation is of some microseconds and the one related to the diffusion of chemical substances (that also represents the material degradation characteristic time itself) is of the order of years. In this case standard incremental techniques must be replaced by other more efficient techniques.

One possibility consists in performing a separated representation of the unknown field:

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{i=N} X_i(\mathbf{x}) \cdot T_i(t) \quad (2)$$

that allows, as we describe later, to non-incremental time integration strategies, which can reduce spectacularly the CPU time.

This space-time separated representation if not a new proposal. In fact such decompositions were proposed many year ago by Pierre Ladeveze as an ingredient of the powerful non-linear-non-incremental LATIN solver that he proposed in the 80s. During the last twenty years many works were successfully accomplished by the Ladeveze's group. The interested reader can refer to [12] [13] [16] and the references therein. In the radial approximation approach (the name given in the pioneer works of Ladeveze) functions depending on space and the ones depending on time were a priori unknown, and they were computed by an appropriate minimization technique.

In this paper we are exploring some new possibilities based on separated space-time representations, that are closely inspired from some existing and well established strategies [4] [12]. In the next section we motivate the use of separated representations and its connection with other more experienced techniques, as the model reduction techniques based on the use of proper orthogonal decompositions. In section 3 we illustrate the application of the Proper Generalized Decomposition through an academic parabolic model. In section 4 we consider the coupling between global and many species kinetic local models and the issue related to the existence of different characteristic times of both the local and global model. This issue was addressed in the context of proper generalized decompositions in [17]. Finally, in section 5 we address a numerical example.

## 2 Motivating the use of separated space-time representations

A well established and widely employed technique allowing to define a separated representation of a given space-time function is based on the application of a proper orthogonal decomposition. We are illustrating the main ideas related to this technique.

Let $u(\mathbf{x}, t)$ be the solution of a certain transient model (in what follows $\mathbf{x} \in \Omega \subset \mathbb{R}^d$, $d = 1, 2, 3$, and $t \in \mathcal{I} \subset \mathbb{R}^+$). We are also assuming that this field is known in a discrete manner, that is, at some points $\mathbf{x}_i$ (the nodes of a
mesh or a grid) and at certain times $t_p = p \times \Delta t$, where $i \in [1, \dots, N_n]$ and $p \in [1, \dots, P]$.

Now, we introduce the notation $u_i^p \equiv u(\mathbf{x}_i, t_p)$ and construct the matrix $\mathbf{Q}$ that contains the snapshots:

$$ \mathbf{Q} = \begin{pmatrix} u_1^1 & u_1^2 & \cdots & u_1^P \\ u_2^1 & u_2^2 & \cdots & u_2^P \\ \vdots & \vdots & \ddots & \vdots \\ u_{N_n}^1 & u_{N_n}^2 & \cdots & u_{N_n}^P \end{pmatrix} \quad (3) $$

The proper orthogonal decomposition (POD) of this discrete field consists in solving the eigenvalue problem:

$$ (\mathbf{Q}\mathbf{Q}^T) \phi = \lambda\phi \quad (4) $$

that results in $N_n$ couples eigenvalue-eigenvector $(\lambda_i, \phi_i)$, $i = 1, \dots, N_n$.

When the field evolves smoothly, the magnitude of the eigenvalues decreases very fast, fact that reveals that the evolution of the field can be approximated from a reduced number of modes (eigenvectors). Thus, if we define a cutoff value $\epsilon$ ($\epsilon = 10^{-8} \times \lambda_1$ in practice, $\lambda_1$ being the highest eigenvalue) only a reduced number of modes are retained. Let $R$ ($R << N_n$) be the number of modes retained, i.e. $\lambda_i \ge 10^{-8} \times \lambda_1$, $i = 1, \dots, R$ and $\lambda_i < 10^{-8} \times \lambda_1, \forall i > R$ (the eigenvalues are assumed ordered). Thus, one could write:

$$ u(\mathbf{x}, t) \approx \sum_{i=1}^{i=R} \phi_i(\mathbf{x}) \cdot T_i(t) \equiv \sum_{i=1}^{i=R} X_i(\mathbf{x}) \cdot T_i(t) \quad (5) $$

where for the sake of clarity the space modes $\phi_i(\mathbf{x})$ will be, from now on, denoted as $X_i(\mathbf{x})$. Eq. (5) represents a natural separated representation (also known as finite sums decomposition).

These modes could be now used to solve other "similar" problems, that is, models involving slight changes in the boundary conditions, model parameters, ... [18] [14] [21]. Other possibility is computing the reduced basis from the standard transient solution within a short time interval (with respect to the whole time interval in which the model is defined) and then solve the remaining part of the time interval by employing the reduced basis. Obviously, both strategies induce the introduction of an error whose evaluation, control and reduction is a challenging issue.

One possibility to construct an adaptive reduced approximation basis, that should be the best reduced approximation basis for the treated problem, consists in alternating a reduction step (based on the application of the proper orthogonal decomposition) and an enrichment stage to improve the quality of the reduced approximation basis in order to capture all the solution features. We proposed recently an enrichment technique based on the use of some Krylov's subspaces generated by the equation residual. This technique known as "a priori" model reduction was originally proposed in [22], widely described in [23]
and successfully applied for solving complex fluid flows within the kinetic theory framework [2] [9] and for speeding up thermomechanical simulations [10]. However, some difficulties were noticed in the application of this strategy: (i) the enrichment based on the use of the Krylov's subspaces is far to be optimal in a variety of models (e.g. the wave equation); (ii) the incremental nature of the algorithm; ...

From the previous analysis we can conclude: (i) the transient solution of numerous models can be expressed using a very reduced number of products each one involving a function of time and a function of space; and (ii) the functions involved in these functional products can be determined simultaneously by applying an appropriate algorithm.

In what follows we are illustrating the simplest strategy able to compute these separated functional couples. As just commented this strategy was proposed by Ladeveze many years ago.

## 3 Illustrating the Proper Generalized Decomposition

In this section we are illustrating the discretization of partial differential equations using a separated representation (radial approximation in the Ladeveze's terminology) of the unknown field.

Let us consider the advection-diffusion equation

$$ \frac{\partial u}{\partial t} - a\Delta u + \mathbf{v} \cdot \nabla u = f(\mathbf{x}, t) \quad \text{in } \Omega \times (0, T_{\text{max}}] \quad (6) $$

with the following initial and boundary conditions

$$ \left\{ \begin{array}{ll} u(\mathbf{x}, 0) = u^0 & \mathbf{x} \in \Omega, \\ u(\mathbf{x}, t) = u_g & (\mathbf{x}, t) \in \partial\Omega \times (0, T_{\text{max}}] \end{array} \right. \quad (7) $$

where $a$ is the diffusion coefficient and $\mathbf{v}$ the velocity field, $\Omega \subset \mathbb{R}^d, d \ge 1$, $T_{\text{max}} > 0$. The aim of the separated representation method is to compute $N$ couples of functions $\{(X_i, T_i)\}_{i=1,...,N}$ such that $\{X_i\}_{i=1,...,N}$ and $\{T_i\}_{i=1,...,N}$ are defined respectively in $\bar{\Omega}$ ($\bar{\Omega}$ being the closure of $\Omega$) and $[0, T_{\text{max}}]$ and the solution $u$ of this problem can be written in the separate form

$$ u(\mathbf{x}, t) \approx \sum_{i=1}^{N} T_i(t) \cdot X_i(\mathbf{x}) \quad (8) $$

The weak form of problem (7) yields:

Find $u(\mathbf{x}, t)$ verifying the boundary conditions (7) such that

$$ \int_{0}^{T_{\text{max}}} \int_{\Omega} u^* \left( \frac{\partial u}{\partial t} - a\Delta u + \mathbf{v} \cdot \nabla u - f(\mathbf{x}, t) \right) d\mathbf{x} dt = 0 \quad (9) $$

for all the functions $u^*(x, t)$ in an appropriate functional space.
We compute now the functions involved in the sum (8). We suppose that the set of functional couples $\{(X_i, T_i)\}_{i=1,...,n}$ with $0 < n < N$ are already known (they have been previously computed) and that at the present iteration we search the enrichment couple $(R(t), S(\mathbf{x}))$ by applying an alternating directions fixed point algorithm that after convergence will constitute the next functional couple $(X_{n+1}, T_{n+1})$. Hence, at the present iteration, $n + 1$, we assume the separated representation

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{n} T_i(t) \cdot X_i(\mathbf{x}) + R(t) \cdot S(\mathbf{x}) \quad (10)$$

The weighting function $u^*$ is then assumed as

$$u^* = S \cdot R^* + R \cdot S^* \quad (11)$$

Introducing (10) and (11) into (9) it results

$$\int_{0}^{T_{max}} \int_{\Omega} (S \cdot R^* + R \cdot S^*) \cdot \left( S \cdot \frac{\partial R}{\partial t} - a \Delta S \cdot R + (\mathbf{v} \cdot \nabla S) \cdot R \right) d\mathbf{x} dt = \\ = \int_{0}^{T_{max}} \int_{\Omega} (S \cdot R^* + R \cdot S^*) \cdot \left( f(\mathbf{x}, t) - \sum_{i=1}^{n} X_i \cdot \frac{\partial T_i}{\partial t} + \right. \\ \left. + a \sum_{i=1}^{n} \Delta X_i \cdot T_i - \sum_{i=1}^{n} (\mathbf{v} \cdot \nabla X_i) \cdot T_i \right) d\mathbf{x} dt \quad (12)$$

We apply an alternating directions fixed point algorithm to compute the couple of functions $(R, S)$:

* Computing the function $S(\mathbf{x})$.

First, we suppose that $R$ is known, implying that $R^*$ vanishes in (11). Thus, Eq. (12) writes

$$\int_{\Omega} S^* \cdot (\alpha_t S - a\beta_t \Delta S + \beta_t \mathbf{v} \cdot \nabla S) d\mathbf{x} = \\ = \int_{\Omega} S^* \cdot \left( \gamma_t(\mathbf{x}) - \sum_{i=1}^{n} \alpha_t^i X_i + a \sum_{i=1}^{n} \beta_t^i \Delta X_i - \sum_{i=1}^{n} \beta_t^i \mathbf{v} \cdot \nabla X_i \right) d\mathbf{x} \quad (13)$$
where

$$ egin{cases} 
\alpha_t = \int_0^{T_{max}} R(t) \cdot \frac{\partial R}{\partial t}(t) dt \\ 
\alpha_t^i = \int_0^{T_{max}} R(t) \cdot \frac{\partial T_i}{\partial t}(t) dt \\ 
\beta_t = \int_0^{T_{max}} R^2(t) dt \\ 
\beta_t^i = \int_0^{T_{max}} R(t) \cdot T_i(t) dt \\ 
\gamma_t(\mathbf{x}) = \int_0^{T_{max}} R(t) \cdot f(\mathbf{x}, t) dt; \quad \forall \mathbf{x} \in \Omega 
\end{cases} 
\quad (14) $$

The weak formulation (13) is satisfied for all $S^*$, therefore we could come back to the associated strong formulation

$$ 
\begin{aligned}
\alpha_t S - a \beta_t \Delta S + \beta_t \mathbf{v} \cdot \nabla S &= \\
&= \gamma_t - \sum_{i=1}^{n} \alpha_t^i X_i + a \sum_{i=1}^{n} \beta_t^i \Delta X_i - \sum_{i=1}^{n} \beta_t^i \mathbf{v} \cdot \nabla X_i
\end{aligned} 
\quad (15) $$

that one could solve by using any appropriate discretization technique for computing the space function $S(\mathbf{x})$.

* Computing the function $R(t)$.

From the function $S(\mathbf{x})$ just computed, we search $R(t)$. In this case $S^*$ vanishes in (11) and Eq. (12) reduces to

$$ 
\begin{aligned}
\int_0^{T_{max}} \int_{\Omega} (S \cdot R^*) \cdot \left( S \cdot \frac{\partial R}{\partial t} - a \Delta S \cdot R + (\mathbf{v} \cdot \nabla S) \cdot R \right) d\mathbf{x} dt &= \\
&= \int_0^{T_{max}} \int_{\Omega} (S \cdot R^*) \cdot \left( f(x,t) - \sum_{i=1}^n X_i \cdot \frac{\partial T_i}{\partial t} + \right. \\
&\qquad \left. + a \sum_{i=1}^n \Delta X_i \cdot T_i - \sum_{i=1}^n (\mathbf{v} \cdot \nabla X_i) \cdot T_i \right) d\mathbf{x} dt
\end{aligned} 
\quad (16) $$

where all the spatial functions can be integrated in $\Omega$. Thus, by using the
following notations

$$
\left\{
\begin{aligned}
\alpha_x &= \int_{\Omega} S(\mathbf{x}) \cdot \Delta S(\mathbf{x}) \, d\mathbf{x} \\
\alpha_x^i &= \int_{\Omega} S(\mathbf{x}) \cdot \Delta X_i(\mathbf{x}) \, d\mathbf{x} \\
\beta_x &= \int_{\Omega} S^2(\mathbf{x}) \, d\mathbf{x} \\
\beta_x^i &= \int_{\Omega} S(\mathbf{x}) \cdot X_i(\mathbf{x}) \, d\mathbf{x} \\
\lambda_x &= \int_{\Omega} S(\mathbf{x}) \cdot (\mathbf{v} \cdot \nabla S(\mathbf{x})) \, d\mathbf{x} \\
\lambda_x^i &= \int_{\Omega} S(\mathbf{x}) \cdot (\mathbf{v} \cdot \nabla X_i(\mathbf{x})) \, d\mathbf{x}, \\
\gamma_x(t) &= \int_{\Omega} S(\mathbf{x}) \cdot f(\mathbf{x}, t) \, d\mathbf{x}; \quad \forall t
\end{aligned}
\right.
\tag{17}
$$

equation (16) reads

$$
\int_0^{T_{max}} R^* \cdot \left( \beta_x \frac{\partial R}{\partial t} + (\lambda_x - a\alpha_x)R - \gamma_x(t) + \sum_{i=1}^n \beta_x^i \frac{\partial T_i}{\partial t} + \sum_{i=1}^n (\lambda_x^i - a\alpha_x^i)T_i \right) dt = 0 \quad (18)
$$

As Eq. (18) holds for all $S^*$, we could come back to the strong formulation

$$
\beta_x \frac{\partial R}{\partial t} = (a \alpha_x - \lambda_x) R + \gamma_x(t) - \sum_{i=1}^{n} \beta_x^i \frac{\partial T_i}{\partial t} + \sum_{i=1}^{n} (a \alpha_x^i - \lambda_x^i) T_i \quad (19)
$$

which is a first order ordinary differential equation that can be solved
easily (even for extremely small times steps) from its initial condition.

These two steps must be repeated until convergence, that is, until verifying
that both functions reach a fixed point. If we denote by $R^{(q)}(t)$ and $R^{(q-1)}(t)$ the
computed functions $R(t)$ at the present and previous iteration respectively, and
the same for the space functions: $S^{(q)}(\mathbf{x})$ and $S^{(q-1)}(\mathbf{x})$, the stoping criterion
used in this work writes:

$$
e = \left\| R^{(q)}(t) \cdot S^{(q)}(\mathbf{x}) - R^{(q-1)}(t) \cdot S^{(q-1)}(\mathbf{x}) \right\|_2 < 10^{-8} \quad (20)
$$

where $10^{-8}$ represents the square root of the machine precision.

This iteration scheme converges very fast (in about 5 iterations), at least for
all the models that we considered in our former works previously referenced.

We denote by $Q_{n+1}$ the number of iterations for solving this non-linear
problem to determine the enrichment couple of functions $X_{n+1}(\mathbf{x})$ and $T_{n+1}(t)$.
After reaching convergence we write $X_{n+1}(\mathbf{x}) = S(\mathbf{x})$ and $T_{n+1}(t) = R(t)$. The enrichment procedure must continue until reaching the convergence of the enrichment global procedure at iteration $N$, when the separated representation of the unknown field writes:

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{N} X_i(\mathbf{x}) \cdot T_i(t) \qquad (21)$$

The more usual global stopping criteria are:

* For models whose exact solution $u^{ref}$ is known:

$$E = \frac{\|u - u^{ref}\|_2}{\|u^{ref}\|_2} < \epsilon \qquad (22)$$

* For models whose exact solution is not known:

$$E = \frac{\left\| \frac{\partial u}{\partial t} - a\Delta u + \mathbf{v} \cdot \nabla u - f(\mathbf{x}, t) \right\|_2}{\|f(\mathbf{x}, t)\|_2} < \epsilon \qquad (23)$$

with $\epsilon$ a small enough parameter ($\epsilon = 10^{-6}$ in our simulations).

**Discussion.** The just proposed strategy needs for the solution of about $N \times Q$ space and time problems (with $Q = (Q_1 + \cdots + Q_N)/N$ and $N$ the number of functional couples needed to approximate, up to the desired precision, the searched solution). Thus, one must compute $N \times Q$ dD problems, $d = 1, 2, 3$, whose complexity depends on the spatial mesh considered, and also $N \times Q$ 1D problems (defined in the time interval $\mathcal{I}$) that only need the solution of an ordinary differential equation from its initial condition. Obviously, even for extremely small time steps, the solution of these transient 1D problems does not introduce any difficulty.

If instead the separated representation just discussed, one performs a standard incremental solution, $P$ dD models, $d = 1, 2, 3$, must be solved ($P$ being the number of time steps, i.e. $P = T_{max}/\Delta t$, where the time step $\Delta t$ must be chosen for verifying the stability conditions).

In all the analyzed cases $N$ and $Q$ are of the order of tens that implies the solution of about hundred three-dimensional problems defined in $\Omega$, instead the thousands (or even millions) needed for solving those models using standard incremental solvers.

A first comparison between both kind of approaches (the one based on the separated representation and the one based on standard incremental strategies) was presented in [5].

In what follows we are introducing new coordinates in the model and by this reason we need to address the solution of more general models that could contain coordinates other than space and time. We proposed in [3] and [4] a technique based on the used of separated representations for treating steady or transient multidimensional models that is a simple generalization of the just presented strategy.
## 4 Efficient coupling of global and local models

We are considering the model given by Eq. (6) in absence of advection, i.e. $v = 0$, and in a one-dimensional physical space $\Omega$:

$$ \frac{\partial u}{\partial t} - a\Delta u = f(x, t) \quad \text{in } \Omega \times (0, T_{\text{max}}] \quad (24) $$

with the following initial and boundary conditions

$$ \left\{ \begin{array}{ll} u(x, 0) = u^0 & x \in \Omega, \\ u(x, t) = u_g & (x, t) \in \partial\Omega \times (0, T_{\text{max}}] \end{array} \right. \quad (25) $$

We are assuming that the term source depends on the local value of $r$ fields $C_i(x, t)$, $i = 1, \dots, r$:

$$ f(x, t) = \sum_{i=1}^{i=r} \gamma_i \cdot C_i(x, t) \quad (26) $$

where the time evolution of the $r$ fields $C_i(x, t)$ is governed by $r$ coupled ordinary differential equations (the so-called kinetic model). For the sake of simplicity we consider the linear case, the non-linear one reduces to a sequence of linear problems by applying an appropriate linearization strategy [5]. The system of linear ODEs writes at each point $x \in \Omega$:

$$ \frac{dC_i(x, t)}{dt} = \sum_{j=1}^{j=r} \alpha_{ij}(x) C_j(x, t) \quad (27) $$

We are assuming that the kinetic coefficients $\alpha_{ij}$ evolves smoothly in $\Omega$, because in practical applications these coefficients depend on the solution of the diffusion problem, $u(x, t)$. For the sake of simplicity and without loss of generality, the alpha coefficients will be assumed later evolving linearly in $x$, but in the description that follows we assume those coefficients constant.

Now, we are describing three possible procedures for solving Eqs. (24) and (27).

1. The simplest strategy consists in using a separated representation of the global problem solution (24) whereas the local problems are integrated in the whole time interval at each nodal position (or integration point). Obviously, this strategy implies the solution of $r$ coupled ordinary differential equations at each node (or integration point). Moreover, the resulting fields $C_i(x, t)$, $r = 1, \dots, r$, don't have a separated structure, and by this reason before to injecting these fields into the global problem (24) we should separate them by invoking, for example, the singular value decomposition (SVD) leading to:

$$ C_i(x, t) \approx \sum_{q=1}^{q=m} X_q^{C,i}(x) \cdot T_q^{C,i}(t) \quad (28) $$
As soon as the source term has a separated structure, the procedure illustrated in the previous sections can be applied again for computing the new trial solution of the global problem that writes:

$$u(\mathbf{x}, t) \approx \sum_{i=1}^{N} X_{i}^{u}(\mathbf{x}) \cdot T_{i}^{u}(t) \qquad (29)$$

Thus, this coupling strategy requires the solution of many local problems (for all the coupled species at all nodal positions (or integration points)). Moreover, after these solutions (that we recall that could be performed in parallel) a singular value decomposition must be applied in order to separate these solutions prior to inject them in the PGD solver of the global problem (24).

2. The second coupling strategy lies in globalizing the solution of the local problems. Thus, we assume that the field related to each species can be written in a separated form:

$$C_i(x, t) \approx \sum_{q=1}^{q=m} X_q^{C,i}(x) \cdot T_q^{C,i}(t) \qquad (30)$$

and now, we apply the procedure described in section 2 to build-up the reduced separated approximation, i.e. for constructing all the functions involved in (30). Thus, instead of solving the $r$ coupled ODEs in Eq. (28) at each nodal position (or integration point), we should solve only $r$ higher dimensional coupled models defined in the physical space and time. Obviously, if the number of nodes (or integration points) is important (mainly when 3D physical spaces are considered) the present coupling strategy could offer significant CPU time savings.

This strategy allows computing directly a separated representation, and then, with respect to the previous one, the application of the SVD is avoided. However, if the number of species is high, the computational efforts can become important, because the space-time separated solver must be applied to each species.

3. The third alternative, that in our opinion is the more appealing one for solving models involving many species, as large as one wants, implies the definition of a new variable $C(x, t, c)$, that as we can notice contains an extra coordinate $c$, with discrete nature, and that takes integer values: $c = 1, \dots, r$, in such manner that $C(x, t, i) \equiv C_i(x, t)$, $i = 1, \dots, r$. Thus, we have increased the dimensionality of the problem, but now, only a single problem should be solved, instead one for each species as was the case when using the previous strategy. This increase of the model dimensionality is not dramatic because as argued in the first section of this work, the separated representation allows circumventing the curse of dimensionality,
allowing for fast and accurate solutions of highly multidimensional models. Now, the issue is the derivation of the governing equation for this new variable $C(x, t, c)$ and the separated representation constructor able to define the approximation:

$$C(x, t, c) \approx \sum_{q=1}^{q=S} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) \quad (31)$$

As this strategy will be retained in our simulations we are focusing in its associated computational aspects in the next section.

## 4.1 Fully globalized local models

The third strategy just referred implies the solution of a single multidimensional model involving the field $C(x, t, c)$. This original introduction deserves some additional comments. The first one concerns the discrete nature of the kinetic equations

$$\frac{dC_i(x, t)}{dt} = \sum_{j=1}^{j=r} \alpha_{ij}(x) \cdot C_j(x, t) \quad (32)$$

Now, by introducing $C(x, t, c)$, such that $C(x, t, i) \equiv C_i(x, t)$, the kinetic equations could be written as:

$$\frac{dC}{dt} = \mathcal{L}_c(C) \quad (33)$$

where $\mathcal{L}_c$ is an operator in the c-coordinate.

If for one second we try to discretize Eq. (33) by finite differences, we could write at each node $(x_k, t_p, i)$:

$$\frac{C(x_k, t_p, i) - C(x_k, t_{p-1}, i)}{\Delta t} = \mathcal{L}_c(C)|_i \quad (34)$$

where

$$\mathcal{L}_c(C)|_i = \sum_{j=1}^{j=r} \alpha_{ij} \cdot C(x_k, t_p, j) \quad (35)$$

represents the discrete form of the c-operator.

Now, we come back to the separated representation constructor for defining the approximation:

$$C(x, t, c) \approx \sum_{q=1}^{q=S} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) \quad (36)$$
For defining such approximation one should repeat the procedure deeply illustrated in section 3. As the operator here involved is less standard we are summarizing the main steps.

We assume that the first $n$ iterations allowed computing the first $n$ sums of Eq. (36)

$$C(x, t, c) \approx \sum_{q=1}^{q=n} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) \quad (37)$$

and now, we look for the enrichment $R(x) \cdot S(t) \cdot W(c)$, such that

$$C(x, t, c) \approx \sum_{q=1}^{q=n} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) + R(x) \cdot S(t) \cdot W(c) = \\ = C^n(x, t, c) + R(x) \cdot S(t) \cdot W(c) \quad (38)$$

satisfies

$$\int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} C^{*}(x, t, c) \cdot \left( \frac{dC}{dt} - \mathcal{L}_{c}(C) \right) dc \, dt \, dx = 0 \quad (39)$$

Obviously, due to the discrete character of the third coordinate, an integration quadrature consisting of $r$ points, $c_1 = 1, \dots, c_r = r$ will be considered later.

Now, for computing the three enrichment functions we are considering again (as in section 3) an alternating directions strategy, that proceeds in three steps (that are repeated until reaching convergence):

1. Assuming functions $S(t)$ and $W(c)$ known, the trial function $C^*(x, t, c)$ writes $R^*(x) \cdot S(t) \cdot W(c)$. Thus the weak form (39) reads:

$$\int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R^* \cdot S \cdot W \cdot (R \cdot S' \cdot W - R \cdot S \cdot \mathcal{L}_c(W)) \, dc \, dt \, dx = \\ - \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R^* \cdot S \cdot W \cdot \sum_{q=1}^{q=n} (X_q^C \cdot (T_q^C)' \cdot A_q) \, dc \, dt \, dx + \\ + \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R^* \cdot S \cdot W \cdot \sum_{q=1}^{q=n} (X_q^C \cdot T_q^C \cdot \mathcal{L}_c(A_q)) \, dc \, dt \, dx \quad (40)$$

where $S' = \frac{dS}{dt}$ and $(T_q^C)' = \frac{dT_q^C}{dt}$.

Now, time integrals and the ones involving the $c$-coordinate can be performed. The ones involving the time are performed as indicated in section 3. The ones involving the coordinate $c$ write:
$$ \int_{0}^{r} W \cdot W \, dc = \sum_{i=1}^{i=r} W(c_i)^2 \quad (41) $$

where as just mentioned $c_i = i, \forall i$,

$$ \int_{0}^{r} W \cdot \mathcal{L}_{c}(W) \, dc = \sum_{i=1}^{i=r} \left( W(c_i) \cdot \sum_{j=1}^{j=r} \alpha_{ij} W(c_j) \right) \quad (42) $$

and similar expressions can be derived for the integrals involved in the right hand member.

Thus, finally it results:

$$ \xi^x \int_{\Omega_x} R^* \cdot R \, dx = \int_{\Omega_x} R^* F^x(x) \, dx \quad (43) $$

where the coefficient $\xi^x$ contains all the integrals in the time and c-coordinates related to the left hand member of Eq. (40) and $F^x(x)$ all the integrals appearing in the right hand member. The strong form related to Eq. (43) writes:

$$ \xi^x R(x) = F^x(x) \quad (44) $$

whose algebraic nature derives from the fact that kinetic model is local and then it does not involve space derivatives.

2. Assuming functions $R(x)$ and $W(c)$ known, the trial function $C^*(x, t, c)$ writes $R(x) \cdot S^*(t) \cdot W(c)$. Thus the weak form (39) reads:

$$ \begin{aligned} & \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S^* \cdot W \cdot (R \cdot S' \cdot W - R \cdot S \cdot \mathcal{L}_c(W)) \, dc \, dt \, dx = \\ & - \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S^* \cdot W \cdot \sum_{q=1}^{q=n} (X_q^C \cdot (T_q^C)' \cdot A_q) \, dc \, dt \, dx + \\ & + \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S^* \cdot W \cdot \sum_{q=1}^{q=n} (X_q^C \cdot T_q^C \cdot \mathcal{L}_c(A_q)) \, dc \, dt \, dx \end{aligned} \quad (45) $$

Now, integrals defined in the physical space $\Omega$ must be computed, but this task does not involve additional difficulties.

Finally it results:

$$ \int_{0}^{T_{max}} S^* \cdot \left( \xi^t S + v^t \frac{dS}{dt} \right) \, dt = \int_{0}^{T_{max}} S^* F^t(t) \, dt \quad (46) $$
where coefficients $\xi^t$ and $v^t$ contain all the integrals in the space and the c-coordinate related to the left hand member of Eq. (45) and $F^t(t)$ the associated integrals appearing in the right hand member. The strong form related to Eq. (46) writes:

$$v^t \frac{dS}{dt} + \xi^t S(t) = F^t(t) \quad (47)$$

whose first order differential nature derives from the form of kinetic models.

3. Assuming functions $R(x)$ and $S(t)$ known, the trial function $C^*(x, t, c)$ writes $R(x) \cdot S(t) \cdot W^*(c)$. Thus the weak form (39) reads:

$$\int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S \cdot W^* \cdot (R \cdot S' \cdot W - R \cdot S \cdot \mathcal{L}_c(W)) \, dc \, dt \, dx = - \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S \cdot W^* \cdot \sum_{q=1}^{q=n} (X_q^C \cdot (T_q^C)' \cdot A_q) \, dc \, dt \, dx + \int_{\Omega} \int_{0}^{T_{max}} \int_{0}^{r} R \cdot S \cdot W^* \cdot \sum_{q=1}^{q=n} (X_q^C \cdot T_q^C \cdot \mathcal{L}_c(A_q)) \, dc \, dt \, dx \quad (48)$$

After performing integration in space and time, it results:

$$\int_{0}^{r} W^{*} \cdot (\xi^{c} W - v^{c} \mathcal{L}_{c}(W)) \, dc = \int_{0}^{r} W^{*} F^{c}(c) \, dc \quad (49)$$

where coefficients $\xi^c$ and $v^c$ contain all the integrals in the space and time related to the left hand member of Eq. (48) and $F^c(c)$ the associated integrals appearing in the right hand member. The strong form related to Eq. (49) writes:

$$-v^c \mathcal{L}_c(W) + \xi^c W(c) = F^c(c) \quad (50)$$

that results in the algebraic system:

$$-v^c \sum_{j=1}^{j=r} \alpha_{ij} W(c_j) + \xi^c W(c_i) = F^c(c_i), \quad i = 1, \dots, r \quad (51)$$
## 4.2 Heterogeneous time integration

In section 3 we illustrated the construction of a separated representation of standard global models, as the one defined in Eq. (24). In this section we illustrated some original possibilities for globalizing local models consisting in kinetic equations. Now, the separated representation of the different kinetic fields $C_i$ can be extracted from the general solution:

$$C(x, t, c) \approx \sum_{q=1}^{q=S} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) \quad (52)$$

by writing:

$$C_i(x, t) = C(x, t, c_i) = C(x, t, i) \approx \sum_{q=1}^{q=S} X_q^C(x) \cdot T_q^C(t) \cdot A_q(i) \quad (53)$$

Thus, we obtain a fast separated representation of the source term in Eq. (24) allowing solving this equation as described in section 3. However, it remains a detail that deserves some comments. The characteristic time related to local and global problems could differ in some orders of magnitude. In that case functions $T_q^C$ are defined using a time step $\delta t$ much lower than the one employed for defining functions $T_i^u$ appearing in the separated representation of the global problem solution $u(x, t)$

$$u(x, t) \approx \sum_{i=1}^{N} X_i^u(x) \cdot T_i^u(t) \quad (54)$$

Thus, the issue is knowing $C_i(x, t)$ with a time resolution of $\delta t$, what is the simplest consistent transfer to the solution of the global problem that is performed with a time resolution of $\Delta t$? It is easy to prove that the simplest choice consistent with a first order discontinuous Galerkin time integration of the ODE involved in the computation of functions $T_i^u(t)$ (that ensures the conservation of the integrated variable) consists of defining $\bar{C}_i(x, t)$ with a resolution of $\Delta t$, such that

$$\bar{C}_i(x, t = n \times \Delta t) = \frac{1}{\Delta t} \int_{(n-1) \times \Delta t}^{n \times \Delta t} C_i(x, t) \, dt \quad (55)$$

In any case, and contrary to the experiences in using standard incremental strategies, the use of the proper generalized decomposition reduced significantly the impact that heterogeneous time integrations have in the total amount of CPU time. This is due to the fact that the time step only affects to the solution of the one dimensional equations solved for computing functions $T_i^u$ and $T_q^C$, and this step is much more faster than the one needed for computing the space functions in the global problem $X_i^u(\mathbf{x})$ that requires the solution of a steady state
2D or 3D problem. Thus, when the differences of the characteristic times of both models is not too large we could perfectly consider the lowest characteristic time step in the integration of both models.

## 4.3 Numerical example

We are illustrating the fully globalization of local kinetic models deeply described in the previous section, because the solution of global models and its coupling do not involves major difficulties, and both aspects were reported in some of our previous works [5] [19].

To enforce a spatial dependance of the kinetic model solution we are considering a kinetic model in which the kinetic coefficients evolves linearly in space. For the sake of simplicity we are also considering a one-dimensional physical space and 10 species (i.e. $r = 10$).

The associated kinetic equations write:

$$ \frac{dC_i(x,t)}{dt} = \sum_{j=1}^{j=r} \alpha_{ij}(x) C_j(x,t) \quad (56) $$

where

$$ \alpha_{ij}(x) = \alpha'_{ij} + x \cdot \alpha''_{ij} \quad (57) $$

The coefficients $\alpha'_{ij}$ and $\alpha''_{ij}$ were defined by

$$ \alpha'_{ij} = \frac{1}{r^2} i \times j - 0.5 \quad (58) $$

$$ \alpha''_{ij} = \frac{1}{r^2} (r + 1 - i) \times (r + 1 - j) - 0.5 \quad (59) $$

The procedure described in section 4 was applied and a separated representation of the field $C(x, t, c)$ was obtained after 11 iterations, i.e. the solution contains 11 sums:

$$ C(x, t, c) \approx \sum_{q=1}^{q=11} X_q^C(x) \cdot T_q^C(t) \cdot A_q(c) \quad (60) $$

Figs. 1-3 depicts the computed functions $X_q^C(x)$, $T_q^C(x)$ and $A_q^C(x)$, that were denoted by F1, F2 and F3 respectively. Obviously, in Fig. 3 only the values of the different curves for $c = 1, 2, \dots, 10$ have sense.

Figures 4 and 5 depict the space-time distribution of species 1 to 5 and 6 to 10 respectively.
<figure><img src="image_1.png" /><figcaption>Figure 1: Separated representation: spatial functions.</figcaption></figure>

<figure><img src="image_2.png" /><figcaption>Figure 2: Separated representation: time functions.</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>Figure 3: Separated representation: functions defined on the species coordinate.</figcaption></figure>

**Discussion.** Because neither the PGD based solver nor the incremental one are computationally optimized, we are not addressing the computing time involved in the simulation. However, we could quantify the advantage in using the PGD based solver by evaluating the number of operations involved in both solutions. We are assuming the simplest numerical choices: explicit incremental integration, ...

We recall the notation previously introduced: $N$ is the number of functional products involved in the separated approximation of field $C(x, t, c)$ being $Q$ the number of iterations of the alternating direction fixed point algorithm needed for computing each one of the functional products. We consider $N_n$ nodes in which the dynamical system must be integrated, being $P$ the number of time steps and $r$ the number of species involved by the kinetic model.

When a standard incremental technique is used (without any kind of computational improvement) the number of operations scales with $N_n \times P \times r$, whereas when using the PGD based solver the number of operations scales with $N \times Q \times (N_n + r^2 + P)$ (this estimation can be easily derived from the algorithm given in section 4.1).

When we consider a simple problem characterized by $N_n = 1000$, $P = 1000$ and $r = 10$, the PGD based solver requires of about $N \times Q \approx 1000$ iterations. Thus, the number of operations involved by the incremental technique is of order $10^7$ whereas the one related to the PGD based solver requires about $10^6$. In this case there is not a clear advantage in using the PGD based solver. Now, if $N_n = 10^6$, $P = 10^6$ and $r = 10$, the incremental solver needs $10^{13}$ operations instead of the $10^9$ involved in the PGD based solver. In this case the CPU time savings is of four orders of magnitude, and this ratio can be increased by assuming more complex simulation scenarios involving more nodes, time steps
<figure><img src="image_4.png" /><figcaption>Figure 4: Space-time distribution of species 1 to 5.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 5: Space-time distribution of species 6 to 10.</figcaption></figure>
or species.

# 5 Conclusions

In this paper we presented some recent developments in the Proper Generalized Decomposition. This technique allows defining non-incremental strategies that can be easily coupled for solving multiscale models. Even if sophisticated and powerful coupling techniques exist for solving non-linear coupled models (the LATIN proposed by Pierre Ladeveze [12] is an excellent candidate) in this paper we applied a simple coupling based in a fixed point that is enough efficient for the models involving parabolic and local kinetics couplings. The main originality of the present work is the globalization of kinetic models as well as the introduction of an additional dimension for accounting for the different species involved by the kinetic model. Thus, instead of solving many local problems (one for each species) at any nodal position (or integration point), we propose a model that only requires the solution of a single problem defined in a multidimensional space that includes the physical space, the time and an extra coordinate related to the different species. The curse of dimensionality is not an issue for Proper Generalized Decompositions because its separated representation nature. In that representations the complexity scales linearly with the dimension of the space instead of the usual exponential growing of standard mesh based strategies. Moreover the decoupling of the time representations allows simple integrators and facilitate the use of different time steps for solving local and global models.

# References

[1] Y. Achdou and O. Pironneau, Computational methods for option pricing, Siam Frontiers in Applied Mathematics, 2005.

[2] A. Ammar, D. Ryckelynck, F. Chinesta and R. Keunings, On the reduction of kinetic theory models related to finitely extensible dumbbells, J. Non-Newtonian Fluid Mech., 134, 2006, 136-147.

[3] A. Ammar, B. Mokdad, F. Chinesta and R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, J. Non-Newtonian Fluid Mech., 139, 2006, 153-176.

[4] A. Ammar, B. Mokdad, F. Chinesta and R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids. Part II: transient simulation using space-time separated representations, J. Non-Newtonian Fluid Mech., 144, 2007, 98-121.
[5] A. Ammar, M. Normandin, F. Daim, D. Gonzalez, E. Cueto, F. Chinesta, Non-incremental strategies based on separated representations: Applications in computational rheology, Communications in Mathematical Sciences, In press.

[6] G. Beylkin and M. Mohlenkamp, Algorithms for numerical analysis in high dimensions, SIAM J. Sci. Com., 26/6, 2005, 2133-2159.

[7] H.J. Bungartz and M. Griebel, Sparse grids, Acta Numerica, 13, 2004, 1-123.

[8] E. Cancès, M. Defranceschi, W. Kutzelnigg, C. Le Bris, Y. Maday, Computational Quantum Chemistry: a primer, Handbook of Numerical Analysis, Vol. X, Elsevier, 2003, pp. 3-270.

[9] F. Chinesta, A. Ammar, A. Falco and M. Laso, On the reduction of stochastic kinetic theory models of complex fluids, Modeling and Simulation in Materials Science and Engineering, 15, 2007, 639-652.

[10] F. Chinesta, A. Ammar, F. Lemarchand, P. Beauchene and F. Boust, Alleviating mesh constraints: Model reduction, parallel time integration and high resolution homogenization, Computer Methods in Applied Mechanics and Engineering, 197/5, 2008, 400-413.

[11] F. Chinesta, A. Ammar and P. Joyot, The nanometric and micrometric scales of the structure and mechanics of materials revisited: An introduction to the challenges of fully deterministic numerical descriptions, International Journal for Multiscale Computational Engineering, In press.

[12] P. Ladeveze, Nonlinear computational structural mechanics, Springer, NY, 1999.

[13] P. Ladeveze, J. Ch. Passieux, D. Neron, The LATIN multiscale computational method and the proper orthogonal decomposition, Computer Methods in Applied Mechanics and Engineering, In press, doi:10.1016/j.cma.2009.06.023

[14] Y. Maday, E. M. Ronquist, The reduced basis element method: application to a thermal fin problem, SIAM J. Sci. Comput., 26/1, 2004, 240-258.

[15] B. Mokdad, E. Pruliere, A. Ammar and F. Chinesta, On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach, Applied Rheology, 17/2, 2007, 1-14.

[16] A. Nouy, P. Ladeveze, Multiscale computational strategy with time and space homogenization: A radial-type approximation technique for solving microproblems, International Journal of Multiscale Computational Engineering, 170/2, 2004.
[17] D. Nron, D. Dureisseix, A computational strategy for thermo-poroelastic structures with a timespace interface coupling, International Journal for Numerical Methods in Engineering, 75/9, 2008, 1053-1084.

[18] H. M. Park and D. H. Cho, The use of the Karhunen-Love decomposition for the modelling of distributed parameter systems, Chem. Engineer. Science, 51, 1996, 81-98.

[19] E. Pruliere, J. Ferec, F. Chinesta, A. Ammar, An efficient reduced simulation of residual stresses in composites forming processes, International Journal of Material Forming, In press.

[20] T. M. Rassias and J. Simsa Finite Sums Decompositions in Mathematical Analysis, John Wiley and Sons Inc., 1995

[21] D. Ryckelynck, L. Hermanns, F. Chinesta and E. Alarcon, An efficient "a priori" model reduction for boundary element models, Engineering Analysis with Boundary Elements, 29, 2005, 796-801.

[22] D. Ryckelynck, A priori hyperreduction method: an adaptive approach, Journal of Computational Physics, 202, 2005, 346-366.

[23] D. Ryckelynck, F. Chinesta, E. Cueto and A. Ammar, On the a priori model reduction: overview and recent developments, Archives of Computational Methods in Engineering, State of the Art Reviews, 13/1, 2006, 91-128.