On the frontier of the simulation world. When models
involve excessive degrees of freedom

Francisco Chinesta, Amine Ammar

► To cite this version:

Francisco Chinesta, Amine Ammar. On the frontier of the simulation world. When models involve excessive degrees of freedom. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2008, 17, pp.583 - 595. 10.3166/remn.17.583-595. hal-01484890

HAL Id: hal-01484890

https://hal.science/hal-01484890

Submitted on 8 Mar 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
On the frontier of the simulation world

When models involve excessive degrees of freedom

Francisco Chinesta* — Amine Ammar**

*Laboratoire de Mécanique des Systèmes et des Procédés – LMSP
UMR 8106 CNRS-ENSAM
151 Boulevard de l’Hôpital, F-75013 Paris
francisco.chinesta@paris.ensam.fr

**Laboratoire de Rhéologie, UMR 5520 CNRS – INPG – UJF
1301 rue de la piscine, BP 53 Domaine universitaire
F-38041 Grenoble cedex 9
Amine.Ammar@ujf-grenoble.fr

ABSTRACT. In the last years, we have assisted to an impressive progression in the numerical modeling capabilities as a result of the progression in computer science but also in the numerical analysis. Thus, new scales have been explored, allowing modeling of richer and finer physical models. In this work we focus on some models encountered in the microscopic description of the physics, all of them with a common particularity: they involve an impressive number of degrees of freedom or are defined in highly multidimensional spaces.

RÉSUMÉ. Ces dernières années, nous avons assisté à une impressionnante progression de la modélisation numérique aussi bien en sciences informatiques qu’en analyse des méthodes numériques. Par conséquent, de nouvelles échelles ont été explorées permettant ainsi la considération de modèles physiques plus riches et plus fins. Dans ce travail, on s’intéresse à certains modèles issus de la description microscopique de la physique. Ces problèmes ont la particularité commune de demander un impressionnant nombre de degrés de liberté ou sont définis dans des espaces de haute dimension.

KEYWORDS: numerical modeling, multi-dimensional models, separated representations, model reduction, nano-physics, statistical mechanics.

MOTS-CLÉS : modélisation numérique, modèles multidimensionnels, séparation de variables, réduction de modèles, nano-science, mécanique statistique.
# 1. Introduction

The fine description of the mechanics and structure of materials at the micro, nano and sub-nanometric scales introduces some specific challenges related to the impressive number of degrees of freedom required or to the highly dimensional spaces in which those models are defined. Despite the fact that spectacular progresses have been accomplished in the context of computational mechanics in the last decade, the efficient treatment of those models, as we describe in the present work, needs further developments.

The brut force approach cannot be considered as a possibility for treating this kind of models. Even some specialists, as the Nobel Prize R.G. Laughlin, affirmed that no computer existing, or that will ever exist, can break the barriers found in quantum mechanics because it is a catastrophe of dimension (Laughlin, 2000). We can understand the catastrophe of dimension by assuming a model defined in a hyper-cube $\Omega$ of dimension $D$: $\Omega = ]-L, L[^D$. Now, if we define a grid to discretize the model, as it is usually performed in the vast majority of numerical methods (finite differences, finite elements, finite volumes, spectral methods etc.), consisting of $N$ nodes on each direction, the total number of nodes will be $N^D$. If we assume that for example $N = 10$ (an extremely coarse description) and $D = 80$ (much lower than the usual dimensions required in quantum or statistical mechanics) the number of nodes in $\Omega$ reaches the astronomical value of $10^{80}$ that represents the presumed number of elementary particles in the universe!

We come back to the analysis of these systems later. In any case, the progress on this field requires further developments on the physical modeling as well as on the introduction of new ideas and methods in the context of computational physics.

In this work we are exploring two model reduction strategies. The first one, based on the use of the Karhunen-Loève decomposition, will be applied for treating models defined in moderate dimensional spaces. The other technique allows treating models defined in highly multidimensional spaces as those encountered in quantum mechanics or in kinetic theory descriptions of simple fluid and gases, complex fluids involving evolving micro or nanostructures, etc.

In *quantum mechanics* one is confronted with the solution of Schrödinger equation giving the electronic and nuclei distribution. If we assume a system composed of $N_p$ particles (electrons and nuclei), the evolution of the joint wavefunction $\Psi = \Psi(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_{N_p}, t)$ is governed by the Schrödinger equation whose dimensionless form in absence of relativistic and spin effects, writes (Cook, 2005):
$$i \frac{\partial \Psi}{\partial t} = - \sum_{p=1}^{p=N_p} \frac{\nabla_p^2 \Psi}{2m_p} + \sum_{p=1}^{p=N_p} \sum_{\substack{k=1 \\ k>p}}^{k=N_p} V_{pk} \Psi \quad [1]$$

where each particle is defined in the whole physical space $\mathbf{x}_j \in \mathbb{R}^3, i = \sqrt{-1}$ and $\hbar$ represents the Planck's constant divided by $2\pi$. The differential operator $\nabla_p^2$ is defined in the conformation space of each particle, i.e. $\nabla_p^2 = \partial^2/\partial x_p^2 + \partial^2/\partial y_p^2 + \partial^2/\partial z_p^2$. The Coulomb's potential accounting for the inter-particles interactions writes:

$$V_{pk} = \frac{q_p \times q_k}{\| \mathbf{x}_p - \mathbf{x}_k \|} \quad [2]$$

where the masses $m_p$ are unity for electrons, the charges $q_j$ are $-1$ for electrons and $+Z_j$ (atomic numbers) for nuclei.

The wavefunction function is then defined in a highly multidimensional space: $\Psi = \Psi(\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_{N_p}, t): \mathbb{R}^{3 \times N_p} \times \mathbb{R}^+ \to \mathbb{C}$. Due to the curse of dimensionality, illustrated above, its direct solution has been only possible for very reduced quantum systems composed of some particles (one or two nuclei and very few electrons). Some approximated techniques have been proposed for solving higher systems (more complex molecules and crystals), being the most widely used the Hartree-Fock (HF) method and the Density Functional Theory (DFT). They made possible the analysis of multi-electronic systems, but nowadays it is mostly accepted that in some particular cases they represent crude approximations. In these cases the interest of solving directly the Schrödinger equation is renewed.

The concern of performing accurate solutions of the Schrödinger equation lies in the possibility to determine quantum-based inter-atomic potentials that could be injected in molecular dynamics simulations, to address some scenarios, badly represented by ordinary pair-wise or N-body inter-atomic potentials, as the ones encountered in the vicinity of interfaces, cracks, defaults...

Other scale that allows continuous descriptions is the one related to statistical mechanics, in which the individuality of each particle is sacrificed in favor of an averaged description that introduces a distribution function giving the probability of finding at a certain point of the space and time, the microstructure described by the so-called conformational coordinates. Thus, for example, when we are considering a suspension of short fibers, the orientation distribution given the fraction of fibers that a certain point $\mathbf{x}$ and time $t$ are oriented on a certain direction $\mathbf{p}$ writes: $\Psi(\mathbf{x}, \mathbf{p}, t): \mathbb{R}^3 \times S(0,1) \times \mathbb{R}^+ \to \mathbb{R}^+$, where $S(0,1)$ represents the surface of the unit sphere spanned by the unit orientation vector $\mathbf{p}$. This model is defined in a
moderate multidimensional space, but other models involve excessive dimensions (as the case of the Schrödinger equation). One example lies in the kinetic theory modeling of macromolecules that are viewed as a Bead-Spring Chain – BSC –. The BSC chain consists of $N+1$ beads connected by $N$ springs. The bead serves as an interaction point with the solvent and the spring contains the local stiffness information depending on local stretching (see Bird et al., 1987, for more details). In this case the microstructure can be described from the distribution function:

$$ \Psi(\mathbf{x}, \mathbf{q}_1, \dots, \mathbf{q}_N, t) : \mathbb{R}^3 \times \mathbb{R}^{3 \times N} \times \mathbb{R}^+ \rightarrow \mathbb{R}^+, \text{ where } \mathbf{q}_i \text{ are the vectors representing the connectors.} $$

In the kinetic theory framework the equation governing the evolution of conformation distribution is known as the Fokker-Planck equation. This equation results from a simple conservation balance of such distribution function:

$$ \frac{\partial \Psi}{\partial t} + \nabla_{\mathbf{x}} \cdot (\mathbf{v} \Psi) = - \sum_{j=1}^{N} \left( \nabla_{\mathbf{q}_j} \cdot (\dot{\mathbf{q}}_j \Psi) \right) \quad [3] $$

In that follows, we are focusing on the numerical issues associated with the numerical solution of the Fokker-Planck equation because the treatment of the Schrödinger equations needs the introduction of additional elements, as the Pauli principle, the Born-Oppenheimer model... whose description is too ambitious for the purpose of the present short paper.

## 2. On the Fokker-Planck numerical solvers

The Fokker-Planck formalism allows describing the evolution of the configuration distribution function, which represents the probability of finding the microstructure in a particular configuration. Obviously this function depends on time, space and configuration coordinates.

The most usual solutions of kinetic theory models associated with complex fluid flows concern the Brownian dynamics approach used in the CONNFESSIT approach (Ottinger, 1992). The stochastic approach was considered in Somasi et al., (2002) for treating MBS (multi-bead-spring) models defined in highly multidimensional conformation spaces. The same idea was considered in the case of short fiber suspensions flows in Chinesta et al., (2003). In that work, the discretization of the advection dominated Fokker-Planck equation governing the fibers orientation was carried out using a particle technique, where the diffusion term was modeled using random walks. The combination of Brownian dynamics with a macroscopic treatment of motion equations is at the basis of the micro-macro approaches deeply reviewed in Keunings (2004). The reader can find in this excellent review numerous and highly valuable descriptions and references.
The main advantage of Brownian dynamics – BD – simulations is that the computational complexity does not scale exponentially (as in the case of mesh-based techniques) with the dimension of the spaces in which the models are defined, because in general, only the moments of the conformation distribution function (and not the distribution function itself) are needed and then computed. In spite of this important advantage, the solutions computed in the framework of BD contain a significant statistical noise (that makes difficult its use in inverse identification or oprimization procedures). This noise can be reduced by increasing the number of realizations of the associated stochastic processes or by using advanced variance reduction techniques.

Thus, when noise-free solutions are required, continuous approximations using either mesh-based or meshless techniques seem to be suitable. Obviously the simulation in the limit case of low or vanishing diffusion effects must be stabilized using a well-established numerical scheme (characteristics based integration, upwinding, etc.). In a pioneering work Chaubal *et al.*, (1997) considered the smooth particle hydrodynamics – a meshless strategy – for solving kinetic theory models, described using the Fokker-Planck formalism, encountered in complex fluid flows modeling. A deterministic particle approach, very close to that proposed in Chaubal *et al.*, (1997), was analyzed in Ammar and Chinesta (2005) and Ammar, Chinesta and Ryckelynck (2006) using smooth particles, but it was noticed that the impact of smoothing on the solution can be significant.

Some attempts of solving the Fokker-Planck equation using a fixed mesh discretisation exist, for example the one performed in Lozinski and Chauviere (2003). The main difficulties of mesh-based approaches are related to the multidimensional character of the problem. Thus, the linear systems obtained after usual implicit or semi-implicit discretizations are extremely large for a practical inversion. On the other hand, explicit discretizations, which do not require matrix inversions, have the constraint of too small time steps. We can affirm that the applicability of standard mesh-based techniques is restricted to low-dimensional conformation spaces.

New advanced numerical strategies for solving accurately and efficiently the Fokker-Planck equation will be described in the next section. The first numerical strategy can be coupled with standard mesh-based discretizations allowing an impressive reduction of the number of degrees of freedom involved by using reduced approximation bases. The reduced approximation basis is constructed by invoking the Karhunen-Loève – KL – decomposition (also known as proper orthogonal decomposition – POD –) in tandem with an approximation basis enrichment algorithm based on the use of some Krylov's subspaces generated by the PDE residual. As this strategy is coupled with standard discretization techniques it can be only used when the kinetic theory model involves a reduced number of dimensions. In highly multidimensional models, standard (mesh-based) discretization techniques fail, and then the just refereed strategy cannot be applied.
When the dimension of the spaces in which the kinetic theory models are described increases significantly, the only possibility of applying deterministic techniques lies in the use of separated representations and tensor product approximation bases. In the next section a new strategy based on these concepts will be described.

# 3. New advanced numerical strategies

## 3.1. Model reduction based on the Karhunen-Loève decomposition

In this section we introduce a numerical approach allowing fast and accurate computations. The idea is very simple: to consider a domain where a certain model is defined and the associated cloud of nodes are able to represent by interpolation the solution everywhere. In general the number of unknowns scales with the number of nodes, and for this reason even if the solution is evolving in time smoothly all the nodes are used for describing it at each time step. In the reduced modeling that we are describing in that follows the numerical algorithm is able to extract the optimal information describing the evolution of the solution in the entire time simulation interval. Thus, the evolution of the solution can be expressed as a linear combination of a reduced number of functions allowing significant CPU time savings.

The extraction of this relevant information is a well known topic based on the application of the proper orthogonal decomposition, also known as Karhunen-Loève decomposition (see Ryckelynck et al., (2006) and the references therein). Usual reduced models perform the simulation of some similar problem or the desired one in a short time interval. From these solutions the Karhunen-Loève decomposition can be applied, that allows to extract the most relevant functions describing the solution evolution (Lorenz, 1956). Now, it is assumed that the solution of a “similar” problem can be expressed using this reduced approximation basis, allowing a significant reduction on the discrete problem size and then significant CPU time savings (Park and Cho, 1996). However, in general the question related to the accuracy of the computed solutions is usually ignored. An original approach combining the model reduction and the control of the solution accuracy was proposed by Ryckelynck (2005), and applied later in different domains (Ryckelynck et al., 2006; Ammar, Ryckelynck et al., 2006).

### 3.1.1. *The Karhunen-Loève decomposition*

We assume known the evolution of a certain field (scalar, vectorial or tensorial) – that from now on and for the sake of simplicity we assume scalar – depending on the physical space **x** and on time *t*, $u(\mathbf{x},t)$. In practical applications, this field is expressed in a discrete form, i.e., it is known at the nodes of a spatial mesh and for some times $u(\mathbf{x}_i, t^p) \equiv u_i^p$. We introduce also the notation:
$u^p(\mathbf{x}) = u(\mathbf{x}, t = p\Delta t) ; \forall p \in [1, \dots, P]$. The main idea of the Karhunen-Loève (KL) decomposition is how to obtain the most typical or characteristic structure $\varphi(\mathbf{x})$ among these $u^p(\mathbf{x})$, $\forall p$. This is equivalent to obtain functions $\varphi(\mathbf{x})$ maximizing $\alpha$:

$$ \alpha = \frac{\sum_{p=1}^{p=P} \left[ \sum_{i=1}^{i=N} \varphi(\mathbf{x}_i) u^p(\mathbf{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\varphi(\mathbf{x}_i))^2} \quad [4] $$

The maximization $\delta\alpha = 0$ leads to the eigenproblem:

$$ \mathbf{k}\varphi = \alpha\varphi \quad [5] $$

where the two points correlation matrix is given by

$$ \mathbf{k}_{ij} = \sum_{p=1}^{p=P} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \Leftrightarrow \mathbf{k} = \sum_{p=1}^{p=P} \mathbf{u}^p(\mathbf{u}^p)^T \quad [6] $$

### 3.1.2. “A posteriori” Reduced modeling

If the evolution of a certain field is known: $u(\mathbf{x}_i, t^p) = u_i^p$, $\forall i \in [1, \dots, N]$, $\forall p \in [1, \dots, P]$, then the matrix $\mathbf{k}$ can be computed and the eigenproblem given by Equation [5] solved. The solution of Equation [5] results in $N$ couples of eigenvalue-eigenvector. However, in a large variety of models involving regular time evolutions of the solution, the magnitude of the eigenvalues decreases very fast, evidencing that the solution evolution can be represented as a linear combination of a reduced number of functions (the eigenvectors related to the highest eigenvalues). In our numerical simulations we consider the eigenvalues ordered $\alpha_1 > \alpha_2 > \dots > \alpha_N$. The $n$ eigenvalues belonging to the interval $\alpha_1 > \dots > \alpha_n$ with $\alpha_n > \alpha_1 \times 10^{-8}$ and $\alpha_{n+1} < \alpha_1 \times 10^{-8}$ are selected, because their associated eigenvectors are enough to represent accurately the entire solution evolution. In a large variety of models $n \ll N$ and moreover $n$ only depends on the regularity of the solution evolution, but neither on the dimension of the physical space (1D, 2D or 3D) nor on the size of the model $N$.

The reduced approximation basis consists of the $n$ eigenvectors $\varphi_1, \dots, \varphi_n$, allowing to define the basis transformation matrix $\mathbf{B}$:

$$ \mathbf{B} = (\varphi_1, \varphi_2, \dots, \varphi_n) \quad [7] $$
whose size is $N \times n$. Thus, the vector containing the unknown field nodal values **u** can be expressed by:

$$ \mathbf{u} = \sum_{i=1}^{n} \phi_i \xi_i(t) = \mathbf{B} \xi(t) \qquad [8] $$

Now, if we consider the linear system of equations resulting from the discretization of a partial differential equation in the form

$$ \mathbf{A} \mathbf{u}^p = \mathbf{f}^{p-1} \qquad [9] $$

where $\mathbf{f}^{p-1}$ accounts for the solution at the previous time step, then, by introducing Equation [8] it results

$$ \mathbf{A} \mathbf{u}^p = \mathbf{f}^{p-1} \Rightarrow \mathbf{A} \mathbf{B} \xi^p = \mathbf{f}^{p-1} \qquad [10] $$

or by multiplying both terms by $\mathbf{B}^T$:

$$ \mathbf{B}^T \mathbf{A} \mathbf{B} \xi^p = \mathbf{B}^T \mathbf{f}^{p-1} \qquad [11] $$

which proves that the final system of equations is of low order, i.e. the dimensions of $\mathbf{B}^T \mathbf{A} \mathbf{B}$ are $n \times n$, with $n \ll N$, and dimensions of $\xi$ and $\mathbf{B}^T \mathbf{f}^{p-1}$ are $n \times 1$.

### 3.1.3. Enriching the reduced approximation basis

The just described strategy allows very fast computation of large size models. However, it is not guaranteed that this reduced basis that has been built from the solution evolving in a short time interval (or from the solution of a “similar model”) remains accurate for describing the solution in the entire simulation interval (or in the actual model). In this manner, if one compute reduced model solutions and wish to keep the confidence on the related solution, an accuracy check would have been performed and an enrichment strategy would have been defined to adapt the reduced approximation basis in order to capture the new events present in the solution evolutions which could not be described accurately from the original reduced approximation basis. For this purpose, Ryckelynck proposed (Ryckelynck, 2005) to start with a low order approximation basis, using some simple functions (e.g. the initial condition in transient problems) or using the eigenvectors of a “similar problem” previously solved or the ones coming from a non-reduced simulation in a short time interval. Then, we compute $S$ time steps of the evolution problem using the reduced model [11] without changing the approximation basis. After these $S$ time steps, the complete discrete system [10] is built-up, and the residual $\mathbf{R}$ evaluated:
$$R = Au^S - f^{S-1} = AB\xi^S - f^{S-1}$$

[12]

If the norm of the residual is small enough, $||R||<\epsilon$, with $\epsilon$ a threshold value small enough, we can continue for other $S$ time steps using the same approximation basis. On the contrary, if the residual norm is too large, $||R||\ge\epsilon$, we must enrich the approximation basis and compute again the last $S$ time steps. This enrichment is built using some Krylov's subspaces, in our simulations the three first subspaces: $B \leftarrow (B, R, AR, A^2R)$. One could expect that the enrichment process is increasing continuously the size of the reduced approximation basis, but in fact, after reaching the convergence, a Karhunen-Loève decomposition is performed on the whole past time interval in order to extract the significant information as well as to define an orthogonal reduced approximation basis. Now, the resulting updated reduced basis is used for computing the next $S$ time steps.

## 3.2. Model reduction based on the use of separated representations

The model reduction technique described in the previous section allows significant computing time savings; however it can be only applied when the associated model can be treated by using standard mesh-based discretizations. Thus, for addressing models defined in highly multidimensional spaces, new reduction strategies whose computational complexity must scale linearly with the dimension of the space should be proposed and checked.

For this purpose we start writing the polynomial approximation of a generic function $T(\mathbf{x}) = T(x_1, x_2, \dots, x_D)$ in the whole domain as:

$$T(\mathbf{x}) \approx \sum_{i=1}^{i=I} \alpha_i T_1^i(x_1) \times \cdots \times T_D^i(x_D) = \sum_{i=1}^{i=I} \alpha_i \prod_{k=1}^{k=D} T_k^i(x_k) \quad [13]$$

The coefficients $\alpha_i$ in the expression [13] are retained in order to consider in the expansion normalized functions $T_k^i(x_k)$. It is also well known that several model solutions can be approximated by a finite, and sometimes so reduced, number of functions products. Expression [13] involves $I \times N \times D$ degrees of freedom instead of the $N^D$ required by the mesh-based discretization techniques.

In that follows we are describing a new advanced technique, in our knowledge never used until now, that combines a separated representation and an adaptation procedure able to build up gradually each product of functions until reaching the convergence. It has some resemblances with the functional approximation used within the LATIN framework (see P. Ladeveze, 1999 and the references therein) as well as with to the ones employed in the post-Hartree-Fock methods (Cancès et al.,
2003). This technique has been successfully applied in a variety of linear, non linear, stationary and non stationary problems inspired from polymer kinetic theory (Ammar, Mokdad et al., 2006; Ammar et al., 2007; Mokdad et al., 2007). For the sake of simplicity we are considering a simple multi-dimensional diffusion problem:

$$ \left\{ \begin{array}{l} \nabla^2 T = f(\mathbf{x}), \mathbf{x}^T = (x_1, \dots, x_D) \in \Omega = ]0, L[ \\ T(\mathbf{x} \in \partial\Omega) = 0 \end{array} \right. \qquad [14] $$

where the general form of the right term is given by

$$ f(\mathbf{x}) \approx \sum_{i=1}^{i=m} a_i F_1^i(x_1) \cdots F_D^i(x_D) = \sum_{i=1}^{i=m} a_i \prod_{k=1}^{k=D} F_k^i(x_k) \qquad [15] $$

expansion that can be performed by using singular value decomposition.

The iteration scheme used to build up the solution proceeds performing a projection and basis enrichment at each iteration:

1. *Projection step.* If we consider the basis function known at the $n^{th}$ iteration:

$$ T(\mathbf{x}) \approx \sum_{i=1}^{i=n} \alpha_i T_1^i(x_1) \cdots T_D^i(x_D) = \sum_{i=1}^{i=n} \alpha_i \prod_{k=1}^{k=D} T_k^i(x_k) \qquad [16] $$

the alpha coefficients can be computed by discretizing the weak formulation of the problem. For this purpose we must define the test functions that in the Galerkin's framework are given by:

$$ T^*(\mathbf{x}) \approx \sum_{i=1}^{i=n} \alpha_i^* T_1^i(x_1) \cdots T_D^i(x_D) = \sum_{i=1}^{i=n} \alpha_i^* \prod_{k=1}^{k=D} T_k^i(x_k) \qquad [17] $$

The resulting linear system has generally a reduced size $n \times n$, reaching its maximum size (just before the convergence) of $I \times I$, that in the case considered results $I \times I \sim 10 \times 10$.

2. *Enrichment step.* Giving the $n$ alpha coefficients, the approximation basis could be enriched by adding a new product of functions $\prod_{k=1}^{k=D} T_k^{n+1}(x_k)$ that needs for the determination of the $D$ involved functions $T_k^{n+1}(x_k)$. For this purpose the trial function:
$$T(\mathbf{x}) = \sum_{i=1}^{i=n} \alpha_i \prod_{k=1}^{k=D} T_k^i(x_k) + \prod_{k=1}^{k=D} R_k(x_k) \quad [18]$$

is injected in the weak formulation, where $R_k(x_k)$ are the unknown quantities of the non linear system obtained, whose size is $D \times N$. The associated test functions are taken, again in the Galerkin's framework, as:

$$T^*(\mathbf{x}) = \sum_{j=1}^{j=D} \left( R_j^*(x_j) \prod_{\substack{k=1 \\ k \neq j}}^{k=D} R_k(x_k) \right) \quad [19]$$

Now, as soon as the functions $R_k(x_k)$ have been determined, the searched functions $T_k^{n+1}(x_k)$ are obtained by normalizing $R_k(x_k)$.

This algorithm has been successfully used to solve models involving hundred dimensions needing of the order of ~$10^{300}$ degrees of freedom if one proceeds in the finite element framework, solution that only needed of around 20 minutes using a standard personal computer!. A multi-dimensional model with known analytical solution was solved to conclude about the solver accuracy. When piecewise linear one-dimensional approximations of all the involved functions were used, a convergence rate of four was noticed, two times higher than the expected rate. In principle the convergence rates could be enhanced by using richer approximations. The separate representation considered in [13] only needs one-dimensional approximations and one-dimensional integrations, because the integral of a product of functions in a hyper-domain can be written as the product of the one-dimensional integrals.

## REMARKS

- if in the resolution of the non-linear model an alternating directions fixed point scheme is used and piecewise linear 1D-interpolations are considered, then the resulting $N \times N$ linear systems to be solved are tri-diagonal;

- the transient models can be solved by assuming that the time is an additional coordinate. However, the stability needs an up-winding of the time derivatives (Ammar et al., 2007);

- the case of non-homogeneous initial or boundary conditions can be addressed by introducing an appropriate change of variables. In some cases the change of variable to be applied is evident, but in others it is a bit more sophisticated. The use of $R$-functions is being analyzed and the first promising results are being obtained,

- some models (as the ones coming from the kinetic theory descriptions of complex fluid flows) are defined in hyper-spaces that involve the physical domain, eventually geometrically complex, and the multidimensional conformational space,
in general geometrically simple. One example concerns the Fokker-Planck equation associated to short fiber suspensions that, as indicated in the first section, makes use of a orientation distribution function $\Psi(\mathbf{x}, \mathbf{p}, t) : \Omega \in \mathbb{R}^3 \times S(0,1) \times \mathbb{R}^+ \rightarrow \mathbb{R}^+.$

The most general and simplest separated representation (that only involves products of functions up to 3D) writes:

$$ \Psi(\mathbf{x}, \mathbf{p}, t) = \sum_{i=1}^{i=\infty} \alpha_i F_i(\mathbf{x}) G_i(\mathbf{p}) J_i(t) \quad [20] $$

## 4. Conclusions

Despite the first promising results that we obtained, and that allow to compute numerical solutions of models until now never solved, there are numerous open questions. Some of them are: (i) stabilization of the advection terms appearing in the Fokker-Planck equation where the usual (non-separable) upwinding techniques fail (the technique based on the use of residual free bubbles seems to be an appealing choice because no knowledge about the advection direction is a priori required); (ii) error estimation and refinement of the 1D descriptions (this could be carried out by using either 1D wavelet approximations in order to take advantage of its multiresolution properties or by invoking again the appealing properties of residual bubbles); (iii) efficient separation of source terms via multi-linear singular value decompositions (iv) multidimensional integration of poorly-separable functions by using a Monte-Carlo technique; (v) use of parallel computing platforms; (vi) use of richer one-dimensional approximations (spectral) to improve – if possible – the convergence rates and (vii) treatment of non-homogeneous boundary and initial conditions...

## 5. References

* Ammar A., Chinesta F., “A particle strategy for solving the Fokker-Planck equation governing the fiber orientation distribution in steady recirculating flows involving short fiber suspensions”, *Lectures Notes on Computational Science and Engineering*, Springer, 43, 2005, p. 1-16.

* Ammar A., Chinesta F., Ryckelynck D., “Deterministic particle approach of Multi-Bead-Spring polymer models”, *European Journal of Computational Mechanics*, vol. 15, n° 5, 2006, p. 481-494.

* Ammar A., Ryckelynck D., Chinesta F., Keunings R., “On the reduction of kinetic theory models related to finitely extensible dumbbells”, *Journal of Non-Newtonian Fluid Mechanics*, 134, 2006, p. 136-147.

* Ammar A., Mokdad B., Chinesta F., Keunings R., “A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids”, *Journal of Non-Newtonian Fluid Mechanics*, 139, 2006, p. 153-176.
Ammar A., Mokdad B., Chinesta F., Keunings R., "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", *Journal of Non-Newtonian Fluid Mechanics*, 144, 2007, p. 98-121.

Bird R.B., Curtiss C.F., Armstrong R.C., Hassager O., *Dynamic of polymeric liquid*, vol. 2, *Kinetic Theory*, John Wiley & Sons, 1987.

Cancès E., Defranceschi M., Kutzelnigg W., Le Bris C., Maday Y., "Computational Quantum Chemistry: a primer", *Handbook of Numerical Analysis*, vol. X, Elsevier, 2003, p. 3-270.

Chaubal C.V., Srinivasan A., Egecioglu O., Leal L.G., "Smoothed particle hydrodynamics techniques for the solution of kinetic theory problems. Part 1: Method", *Journal of Non-Newtonian Fluid Mechanics*, 70, 1997, p. 125-154.

Chinesta F., Chaidron G., Poitou A., "On the solution of the Fokker-Planck equations in steady recirculating flows involving short fiber suspensions", *Journal of Non-Newtonian Fluid Mechanics*, 113/2-3, 2003, p. 97-125.

Cook D.B., *Handook of computational quantum chemistry*, Dover, 2005.

Keunings R., "Micro-macro methods for the multiscale simulation of viscoelastic flow using molecular models of kinetic theory", *Rheology Reviews 2004*, D.M. Binding and K. Walters (Eds.), British Society of Rheology, 2004, p. 67-98.

Ladeveze P., *Nonlinear computational structural mechanics*, Springer, NY, 1999.

Laughlin R.B., "The theory of everything", *Proceedings of the USA National Academy of Sciences*, 97, 2000.

Lorenz E.N., Empirical orthogonal functions and statistical weather prediction, MIT, Department of Meteorology, Scientific Report N1, Statistical Forecasting Project, 1956.

Lozinski A., Chauviere C., "A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model", *J. Computational Physics*, 189, 2003, p. 607-625.

Mokdad B., Pruliere E., Ammar A., Chinesta F., "On the simulation of kinetic theory models of complex fluids using the Fokker-Planck approach", *Applied Rheology*, vol. 17, n° 2, 26494, 2007, p. 1-14.

Ottinger H.C., *Stochastic processes in polymeric fluids*, Springer, 1996.

Park H.M, Cho D.H., "The use of the Karhunen-Loève decomposition for the modelling of distributed parameter systems", *Chem. Engineer. Science*, 51, 1996, p. 81-98.

Ryckelynck D., "A priori hyperreduction method: an adaptive approach", *Journal of Computational Physics*, 202, 2005, p. 346-366.

Ryckelynck D., Chinesta F., Cueto E., Ammar A., "On the a priori model reduction: overview and recent developments", *Archives of Computational Methods in Engineering*, State of the Art Reviews, vol. 13, n° 1, 2006, p. 91-128.

Somasi M., Khomami B., Woo N.J., Hur J.S., Shaqfeh E.S.G., "Brownian dynamics simulation of bead-rod and bead-spring chains: numerical algorithms and coarse graining issues", *Journal of Non-Newtonian Fluid Mechanics*, 108, 2002, p. 227-255.