# Coupling Finite Elements and Reduced Approximation Bases

Amine Ammar, Etienne Prulière, Julien Férec, Francisco Chinesta, Elías Cueto

► To cite this version:

Amine Ammar, Etienne Prulière, Julien Férec, Francisco Chinesta, Elías Cueto. Coupling Finite Elements and Reduced Approximation Bases. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2009, 18 (5-6), pp.445-463. 10.3166/ejcm.18.445-463. hal-01006702

HAL Id: hal-01006702

https://hal.science/hal-01006702

Submitted on 6 Jan 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Coupling finite elements and reduced approximation bases

Amine Ammar* — Etienne Pruliere** — Julien Férec**
Francisco Chinesta** — Elias Cueto***

*Laboratoire de Rhéologie: UMR CNRS-UJF-INPG
13 rue de la Piscine, BP 53 Domaine Universitaire
F-38041 Grenoble cedex 9, France
Amine.Ammar@ujf-grenoble.fr

**EADS Corporate Foundation International Chair
GeM, UMR CNRS-Ecole Centrale de Nantes, France
1 rue de la Noe, BP 92101, F-44321 Nantes cedex 3, France
Francisco.Chinesta@ec-nantes.fr

***I3A, Universidad de Zaragoza
María de Luna, 7, E-50018 Zaragoza, Spain
ecueto@unizar.es

**ABSTRACT.** Models encountered in computational physics and engineering, usually involve too many degrees of freedom, too many simulation time-steps, too many iterations (e.g. non-linear models, optimization or inverse identification...), or simply excessive simulation time (for example when simulation in real time is envisaged). In some of our former works different reduction techniques were developed, some of them based on the use of an adaptive proper orthogonal decomposition and the other ones based on the use of separated representations. In this paper we are analyzing the coupling between reduced basis and standard finite element descriptions.

**RÉSUMÉ.** Les modèles que l'on retrouve en physique et ingénierie font souvent appel à un nombre excessif de degrés de liberté, de pas de temps, d'itérations (par exemple dans le cas de la résolution de modèles non linéaires, dans l'optimisation ou encore dans l'identification inverse), ou tout simplement d'un temps de calcul excessif (par exemple quand l'on s'intéresse à la simulation en temps réel). Dans certains de nos travaux récents nous avons proposé l'emploi de bases réduites, construites à partir de l'emploi de la décomposition orthogonale aux valeurs propres (POD) ou d'une représentation séparée. Dans ce papier nous nous intéressons au couplage des bases réduites avec des approximations standard de type éléments finis.

**KEYWORDS:** model reduction, proper orthogonal decomposition, separated representations, finite elements.

**MOTS-CLÉS :** réduction de modèles, POD, séparation de variables, éléments finis.
# 1. Introduction

Most engineering systems can be modeled by a continuous model usually expressed as a system of linear or non-linear coupled partial differential equations describing the different conservation balances (momentum, energy, mass and chemically reacting substances). From a practical point of view, the determination of the problem solution, that is, the knowledge of the different fields characterizing the physical system at any point and time (velocity, pressure, temperature, concentrations...) is not possible in real systems due to the complexity of models, geometries and boundary conditions. For this reason the solution is searched only at some points and at some times, from which it could be interpolated to any other point and time. Techniques leading to this kind of representation are known as discretization techniques. There are today numerous discretization techniques (finite elements, finite volumes, finite differences, meshless techniques...). The optimal technique to be applied depends on the model and on the domain geometry. The development of numerical analysis and the computation availabilities make possible today the resolution of complex systems involving millions of unknowns related to the discrete model. However, the complexity of the models is also increasing exponentially, and today engineers are not only interested in solving models, but also in solving these modes very fast and, if possible, with a great accuracy.

In the context of control, optimization, inverse analysis, and in general computation in real time, it is clear that numerous problems must be solved, and for this reason the question related to the computation time becomes crucial. The question is very simple: is it possible to perform very fast and accurate simulations? Different answers have been given to this question depending on the scientific community to which this question is addressed. For specialists in computational science an answer to this question requires the improvement of computational resources, high performance computing and the use of parallel computing platforms. For some specialists in numerical analysis the challenge is in the fast resolution of linear systems via the use of preconditioners or multigrid techniques among many others. For others the idea is to adapt the cloud of nodes (points where the solution is computed) in order to avoid excessive number of unknowns. However, today all these approaches allow to alleviate the computation efforts but the fast and accurate computation remains a real challenge.

This paper presents another different approach, allowing impressive reduction in the number of degrees of freedom without compromising the computed solution accuracy, based on the use of reduced approximation bases. The idea is very simple. Consider a domain where a certain model is defined and the associated cloud of nodes 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 this paper 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 (defining the reduced approximation basis), and then the size of the resulting linear problems are very low, and the CPU time savings of some orders of magnitude can be attained (Ammar *et al.*, 2006).

The extraction of this relevant information is a well known topic that can be
performed by applying the proper orthogonal decomposition – POD –, also known
as Karhunen-Loève decomposition – KLD – (Karhunen, 1946; Loève, 1963) that
will be summarized in the next section. This kind of approach has been widely used
for weather forecast purposes (Lorenz, 1956), turbulence (Sirovich, 1987; Holmes *et
al.*, 1997), solid mechanics (Krysl *et al.*, 2001) but also in the context of chemical
engineering for control purposes (see the Park's works, e.g. Park and Cho, 1996).

Usual reduced models perform the simulation of some similar problems or the
desired one in a short time interval. From these solutions the POD decomposition can
be performed, allowing extracting the most relevant functions describing the solution
evolution. Now, it is assumed that the solution of a “similar” problem can be expressed
using this reduced approximation basis, allowing to a significant reduction on the
discrete problem size and then to significant CPU time savings. 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; Chinesta *et al.*, 2008).

In the adaptive technique initially proposed by Ryckelynck (2005), the so called
“a priori” model reduction, the solution of a transient model is expressed, under
certain regularity requirements, from a finite sum decomposition, that is, as the sum
of some products, each one involving a function of space and a function of time. For
this purpose, one must define a procedure to define the functions of space, and the
ones depending on time are computed by enforcing the verification of the dynamical
system. The main issue of this kind of technique is precisely the procedure to build-
up the optimal reduced approximation basis. Ryckelynck proposed the use of the
residual or some Krylov’s subspaces generated by the residual to enrich the
approximation basis, and after convergence, the just enriched basis is reduced by
applying the proper orthogonal decomposition (the algorithm will be summarized
later). However, some times (e.g. the wave equation with a moving load) the
enrichment procedure based on the use of the Krylov’s subspaces doesn’t run.

Another possibility to define a separated representation of a solution consists of
building-up each term in the finite sum decomposition from the Galerkin form of the
governing equation by using an alternating direction scheme, looking for the fixed
point of the iteration solver. This strategy was originally introduced by Ladeveze
(see Ladeveze 1999 and the references therein) that called it “radial approximation”.
This strategy was later and independently used for the solution of stochastic models
(Nouy, 2009) and by Ammar, Mokdad et al. (2006 and 2007) for addressing highly
multidimensional models. The main advantage of this technique with respect to the
one based on the use of the adaptive POD is that both functions (the one defined in space and the one defined on time) are built-up by the algorithm without any major assumption in a transparent way for the user. Moreover, the first numerical experiments and the first numerical analysis results proved its superior robustness.

Thus, both techniques just presented can be used for defining reduced approximation bases able to compute fast and accurately the solution of transient models. However, some times the solutions of general models could involve localization, singularities, discontinuities (fixed or evolving)... and then the coupling between reduced approximation bases (for describing the smooth part of the solution) and standard finite elements for describing the localized behavior seems a very appealing strategy in computational mechanics. This work is a first attempt towards this ambitious objective.

Even if significant computing time savings are expected when using this coupled strategy, this paper only focus in the analysis of the coupling strategy and the evaluation of the reduction of the number of degrees of freedom involved in the resulting discrete model. A deep analysis of the computing time savings will be addressed in an oncoming paper.

## 2. A priori model reduction based on the use of the POD

We assume that the evolution of a certain field depending on the physical space $\mathbf{x}$ and on time $t$, $u(\mathbf{x},t)$ is known. In practical applications, this field is expressed in a discrete form, that is, it is known at the $N$ nodes of a spatial mesh at $P$ different times $u(\mathbf{x}_i, t^p) \equiv u_i^p$. We can also write introducing a spatial interpolation $u^p(\mathbf{x}) \equiv u(\mathbf{x}, t = p\Delta t)$; $\forall p \in [1, \dots, P]$. The main idea of the POD 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 obtaining functions $\varphi(\mathbf{x})$ maximizing $\alpha$

$$ \alpha = \frac{\sum_{p=1}^{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 [1] $$

which leads to:

$$ \sum_{p=1}^{P} \left[ \left[ \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) u^p(\mathbf{x}_i) \right] \left[ \sum_{j=1}^{j=N} \varphi(\mathbf{x}_j) u^p(\mathbf{x}_j) \right] \right] = \alpha \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) \varphi(\mathbf{x}_i) \quad [2] $$
$\forall \tilde{\phi}$, which can be rewritten in the form

$$ \sum_{i=1}^{i=N} \left[ \sum_{j=1}^{j=N} \left\{ \sum_{p=1}^{p=P} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \varphi(\mathbf{x}_j) \right\} \tilde{\varphi}(\mathbf{x}_i) \right] = \alpha \sum_{i=1}^{i=N} \tilde{\varphi}(\mathbf{x}_i) \varphi(\mathbf{x}_i) \quad [3] $$

Defining the vectors **a** such that its *i*-component is $a(\mathbf{x}_i)$, Equation [3] takes the following matrix form

$$ \tilde{\phi}^T \mathbf{k} \varphi = \alpha \tilde{\phi}^T \varphi ; \quad \forall \tilde{\phi} \quad \Rightarrow \quad \mathbf{k} \varphi = \alpha \varphi \quad [4] $$

where the two points correlation matrix **k** 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 [5] $$

which is symmetric and positive definite. If we define the matrix **Q** containing the discrete field history:

$$ \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^1 & u_N^2 & \cdots & u_N^P \end{pmatrix} \quad [6] $$

is easy to verify that the matrix **k** in Equation [5] results

$$ \mathbf{k} = \mathbf{Q} \cdot \mathbf{Q}^T \quad [7] $$

## 2.1. Reduced modeling

If the evolution of a certain field is known

$$ u(\mathbf{x}_i, t^p) \equiv u_i^p, \quad \forall i \in [1, \dots, N], \quad \forall p \in [1, \dots, P] \quad [8] $$

which is coming from some direct simulations or from experimental measures, then matrices **Q** and **k** can be computed and the eigenvalue problem given by Equation [4] solved. The solution of Equation [4] results in *N* couples of eigenvalue-eigenvector. However, in a large number of models involving smooth 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 ones related to the highest eigenvalues).

In our numerical applications 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 \cdot 10^{-8}$ and $\alpha_{n+1} < \alpha_1 \cdot 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, and neither on the dimension of the physical space (1D, 2D or 3D) nor on the size of the model $N$. We verified in the numerical examples addressed in this work but also in many other addressed in former works that by considering a threshold value of $10^{-8}$ the error between the reduced solution and fully finite element solution ($L^2$ norm) was in the order of $10^{-4}$. If the threshold is reduced to $10^{-2}$ for example, error of about 0.1 ($L^2$ norm) were noticed in the examples addressed in this paper. Even if we cannot conclude, some numerical experiments seem to indicate that the error scales with the square root of the threshold value.

The reduced approximation basis consists on the $n$ eigenvectors $\phi_1, \cdots, \phi_n$, allowing to define the basis transformation matrix $\mathbf{B}$:

$$ \mathbf{B} = (\phi_1, \phi_2, \cdots, \phi_n) \qquad [9] $$

whose size is $N \times n$. Thus, the vector containing the field nodal values $\mathbf{u}$ can be expressed by:

$$ \mathbf{u} = \sum_{i=1}^{n} \phi_i \quad a_i(t) = \mathbf{B} \mathbf{a}(t) \qquad [10] $$

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

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

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

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

that premultiplying both terms by $\mathbf{B}^T$ writes

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

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 the dimensions of both $\mathbf{a}$ and $\mathbf{B}^T \mathbf{f}^{p-1}$ are $n \times 1$.

**REMARK 1.** Equation [13] can be also derived introducing the approximation [10] into the partial differential equation Galerkin form.

## 2.2. Enriching the reduced approximation basis

The just described strategy allows for very fast computation of large size models. For example one could solve the full model using some standard discretization
technique (finite differences, finite elements...) for a small time interval and then define matrix **Q** and **k** allowing compute the reduced approximation basis transformation **B** that leads to the reduced solution procedure illustrated by Equation [13].

However, it is not guaranteed that this reduced basis that has been built from the solution evolving in a short time interval remains accurate for describing the solution in the entire simulation interval. Moreover, during the simulation material properties, boundary conditions... could change, and sometimes these changes are significant. Thus, in general the level of confidence expected for this reduced model solution is decreasing as the time (out the interval that served to build the reduced approximation basis) increases.

In this manner, if one would compute reduced model solutions and keep the confidence on the computed solution, a check of accuracy must be applied and, if needed, an enrichment strategy must be performed in order to adapt the reduced approximation basis for capturing the new events present in the solution evolution which can not be described accurately from the original reduced approximation basis. Using this strategy one can attain computing time savings of some orders of magnitude (Ammar *et al.*, 2006).

For this purpose, Ryckelynck (2005) proposed 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 full simulation within a short time interval. Now, we compute *S* iterations of the evolution problem using the reduced model [13] without changing the approximation basis. After these *S* iterations, the complete discrete system [12] is constructed, and the residual **R** evaluated:

$$ \mathbf{R} = \mathbf{A}\mathbf{u}^S - \mathbf{f}^{S-1} = \mathbf{A}\mathbf{B}\mathbf{a}^S - \mathbf{f}^{S-1} \quad [14] $$

If the norm of the residual is small enough, $||\mathbf{R}|| < \varepsilon$, with $\varepsilon$ a threshold value small enough, we can continue for other *S* iterations using the same approximation basis. On the contrary, if the residual norm is too large, $||\mathbf{R}|| \ge \varepsilon$, we need to enrich the approximation basis and compute again the last *S* iterations. This enrichment is built using some Krylov's subspaces, in our case the three first subspaces: $\mathbf{B} \leftarrow (\mathbf{B}, \mathbf{R}, \mathbf{AR}, \mathbf{A}^2\mathbf{R})$.

In the numerical example here addressed we considered 3 Krylov's subspaces (this choice is not essential) and start by putting $S=1000$. If an enrichment is needed after the *S* iterations, so we put $S=S/2$. If it is not the case, we put: $S=2S$.

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 POD decomposition is performed on the whole past time interval in order to extract the significant information. The reduced basis is then updated accordingly.
Obviously, because this control and enrichment stages we can ensure the quality of the computed solution.

We described in this section the main ideas of the adaptation strategy that we are using. However, in order to be efficient from a computational point of view, the numerical algorithm that we are using is a bit more complex. Interested readers can refer to Ryckelynck et al. (2006) and the references therein for a more detailed and valuable description.

# 3. Towards a finite elements-reduced approximation bases coupling

## 3.1. Generic form of the model

In what follows we assume the physical model described by a generic partial differential equation

$$ \mathcal{K}(u(\mathbf{x},t)) + \mathcal{L}(u(\mathbf{x},t)) = 0 \quad \text{in } [\Omega \subset \mathbb{R}^d] \times [I \subset \mathbb{R}^+] \quad [15] $$

where $d = 1, 2, 3$ and $\mathcal{K}$ and $\mathcal{L}$ are two differential operators. In what follows, and without loss in generality, we are assuming that $\mathcal{K} = \frac{\partial}{\partial t}$. We introduce the following notation:

$$ \left\{ \begin{array}{l} \mathbb{N} = \int_{\Omega} \mathbb{N} \mathbb{N}^T d\Omega \\ \mathbb{L} = \int_{\Omega} \mathbb{N} \mathcal{L} \mathbb{N}^T d\Omega \end{array} \right. \quad [16] $$

where $\mathbb{N}$ is the vector containing the shape functions and matrix $\mathbb{L}$ incorporates integration by parts and stabilizations of the advective terms eventually involved by the operator $\mathcal{L}$.

With this notation the discrete model results:

$$ \mathbf{u}^{*T} \mathbb{N} \dot{\mathbf{u}} + \mathbf{u}^{*T} \mathbb{L} \mathbf{u} = 0, \quad \forall \mathbf{u}^* \Rightarrow \mathbb{N} \dot{\mathbf{u}} + \mathbb{L} \mathbf{u} = \mathbf{0} \quad [17] $$

## 3.2. Defining the approximation coupling

This section deals with models involving localized singularities or discontinuities, being the solution far from these regions quite smooth. When these regions evolve in time, additional difficulties appear whose treatment is out of the scope of the present paper. Thus, one expects that efficient simulations could be achieved by coupling reduced descriptions in the smooth regions and a fully finite element description within the regions involving singularities or discontinuities.
Discontinuities can be efficiently captured in the extended finite elements framework, however the application of the techniques previously described in the context of the extended finite elements for describing evolving discontinuities introduces serious difficulties because the number of approximation functions needed becomes too large. For these reasons the simplest and most efficient solution consists of dividing the domain in two parts, and then applying different approximations in each of these resulting regions. The one that concerns the finite element approximation must cover the domain in which the singularity or discontinuity is evolving. Obviously, if the singularity (or discontinuity) explores the whole domain the technique that we are describing doesn't work, but in the most of engineering models these singularities remain confined in regions in general much smaller than the whole domain. Another difficulty related to the application of the coupling that we are describing is the necessity of an “a priori” knowledge about the region in which the evolving singularity (or discontinuity) will be confined.

Despite all the criticisms that we just addressed to the coupling strategy that we are describing in this section, numerous models accept this simple coupling, as we illustrate in the next section.

Let $\Omega_2$ be the small region of interest where localized singularities exist and $\Omega_1$ the remaining part of the domain (where the solution is assumed varying smoothly):

$$ \begin{cases} \Omega_1 \cup \Omega_2 = \Omega \\ \Omega_1 \cap \Omega_2 = \emptyset \end{cases} \qquad [18] $$

Each node of the discretization belongs to a single subdomain. The matrix form of the discrete model [17] is now rewritten, by introducing a decomposition of the unknown vector $\mathbf{u}$ in $\mathbf{u}_1$ and $\mathbf{u}_2$, containing the degrees of freedom in the domains $\Omega_1$ and $\Omega_2$ respectively, as:

$$ \begin{bmatrix} \mathbb{N}_{11} & \mathbb{N}_{12} \\ \mathbb{N}_{21} & \mathbb{N}_{22} \end{bmatrix} \begin{pmatrix} \dot{\mathbf{u}}_1 \\ \dot{\mathbf{u}}_2 \end{pmatrix} + \begin{bmatrix} \mathbb{L}_{11} & \mathbb{L}_{12} \\ \mathbb{L}_{21} & \mathbb{L}_{22} \end{bmatrix} \begin{pmatrix} \mathbf{u}_1 \\ \mathbf{u}_2 \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ \mathbf{0} \end{pmatrix} \qquad [19] $$

On the domain $\Omega_1$ the “a priori” model reduction described in the previous section could be applied. A set of reduced basis function are involved in the matrix $\mathbf{B}$. For this purpose we introduce the following notation:

$$ \left\{ \begin{aligned} \mathbf{n}_{11} &= \mathbf{B}^T \mathbb{N}_{11} \mathbf{B} \\ \mathbf{l}_{11} &= \mathbf{B}^T \mathbb{L}_{11} \mathbf{B} \\ \mathbf{u}_1 &= \mathbf{B} \mathbf{a}_1 \end{aligned} \right. \qquad [20] $$

which allows writing the system [19] as:
$$ \left( \begin{matrix} \mathbf{n}_{11} & \mathbf{B}^T \mathbf{N}_{12} \\ \mathbf{N}_{21} \mathbf{B} & \mathbf{N}_{22} \end{matrix} \right) \left( \begin{matrix} \dot{\mathbf{a}}_1 \\ \dot{\mathbf{u}}_2 \end{matrix} \right) + \left( \begin{matrix} \mathbf{l}_{11} & \mathbf{B}^T \mathbf{L}_{12} \\ \mathbf{L}_{21} \mathbf{B} & \mathbf{L}_{22} \end{matrix} \right) \left( \begin{matrix} \mathbf{a}_1 \\ \mathbf{u}_2 \end{matrix} \right) = \left( \begin{matrix} \mathbf{0} \\ \mathbf{0} \end{matrix} \right) \quad [21] $$

Thus, at least in the region $\Omega_1$ the computing cost has been significantly reduced. The residual calculation is performed over the whole domain $\Omega$ and for the enrichment of the basis, only the values of the residual on $\Omega_1$ are retained.

Because the use of a finite element interpolation of both the unknown function in region $\Omega_2$ (fully finite element description) and of the functions involved in the reduced approximation basis, the approximation continuity in the whole domain is ensured as when one uses a fully finite element description in the entire domain.

From an algebraic point of view, the band size increases as can be noticed in Equation [21]. However, as the reduced model contains very few degrees of freedom its impact on the linear solver is negligible.

# 4. Numerical examples

## 4.1. A simple heat transfer problem

To illustrate the just described technique we first consider a simple 2D heat transfer problem. This problem concerns the evolution of the temperature field $T(\mathbf{x}, t)$ governed by the heat equation:

$$ \frac{\partial T}{\partial t} - \alpha \Delta T = h(\mathbf{x}), \quad (\mathbf{x}, t) \in \Omega \times I \quad [22] $$

where the diffusivity $\alpha$ is assumed isotropic and constant in the whole domain and time (it was fixed to a unit value in the simulations here reported), $\Omega = (-\pi/2, \pi/2)$ and $I = (0,1]$.

The heat source $h(\mathbf{x})$ was taken in the form

$$ h(x, y) = \max(x^9, 0) \cdot \max(y^9, 0) + \min(x^9, 0) \cdot \min(y^9, 0) \quad [23] $$

in order to induce high gradient of the temperature field at the boundaries neighborhood.

Finally, the initial and boundary conditions are given by

$$ \begin{cases} T(x, y, t = 0) = \cos^2(x) \cos^2(y) \\ T(\mathbf{x} \in \partial \Omega, t) = 0 \end{cases} \quad [24] $$
This problem was solved by employing a reduced approximation in the whole domain and also by applying the coupling technique described in the previous section. Figure 1 illustrates the solution computed by using the first strategy. The solution was approximated with only six couples of functions according to:

$$T(\mathbf{x}, t) = \sum_{i=1}^{i=6} F_i(\mathbf{x}) a_i(t) \qquad [25]$$

whose discrete counterpart writes:

$$\mathbf{T}(t) = \sum_{i=1}^{i=6} \mathbf{F}_i a_i(t) \qquad [26]$$

Figure 1 depicts the four most representative modes ($F_1(\mathbf{x}), \cdots, F_4(\mathbf{x})$), the evolution of all the reduced coordinates $a_i(t)$, $i = 1, \cdots, 6$ as well as the resulting temperature field at $t = 1$. There are no appreciable differences between the finite element solution and the reduced one (the norm of the error is $2 \cdot 10^{-4}$).

The solution employing the coupling method was performed by decomposing the physical domain into two regions, as depicted in Figure 2. In that figure the triangular elements located in the white region have nodes taking part of the standard finite element approximation and others that are concerned by the reduced approximation.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Solution computed using a fully reduced approximation basis: (a)… (d) most representative modes; (e) reduced coordinates evolution; (f) solution at $t = 1$</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** (left) Mesh considered for the functional representation of both the reduced approximation functions and the finite element interpolation; (right) Domain decomposition: the green region consists in all the nodes where the solution is approximated employing a reduced basis whereas the region in red denotes the region in which a fully finite element approximation is used</figcaption></figure>

The solution in the green region is performed by applying the model reduction previously described. To represent this solution only six modes were needed. The four most representative modes, the evolution of the reduced coordinates and the solution at time $t = 1$ are depicted in Figure 3. Finally Figure 4 depicts the finite element solution at six different times. By representing together the solutions at any time computed in the green and red regions, the solution at that time is defined in the whole domain.

## 4.2. Accounting for localized behaviour in kinetic theory models

The description of the microstructural conformation in complex fluids needs too many degrees of freedom. In this kind of models one must describe the molecular conformation that in the simplest descriptions (that correspond to the most widely used) is carried out by defining the molecule end-to-end vector **x**. Thus, a molecule is finally modelled as two beads connected by a spring. When this spring is Hookean it allows infinite extensions, however in order to limit the maximum molecule stretching different non-linear spring stiffness were introduced. One of the most popular is the so-called FENE model, in which the maximum extension is denoted by $\sqrt{b}$ and the spring stiffness writes:

$$f(\mathbf{x}) = \frac{1}{1 - \frac{\mathbf{x}^2}{b}} \qquad [27]$$
The microstructure state is then defined by introducing a probability distribution function (p.d.f.) $\Psi(\mathbf{x}, t)$ giving the fraction of molecules that at time $t$ has the end-to-end vector given by $\mathbf{x}$. Obviously, because the finite extension of such molecules $\mathbf{x} \in B(0, \sqrt{b})$, that is, the end-to-end vector belongs to the ball of radius $\sqrt{b}$.

The dimensionless conservation balance of the p.d.f. $\Psi(\mathbf{x}, t)$ writes

$$ \mathcal{K}(\Psi(\mathbf{x}, t)) + \mathcal{L}(\Psi(\mathbf{x}, t)) = 0 \quad \text{in} \quad [\Omega = B(0, \sqrt{b})] \times [I \subset \mathbb{R}^+] \quad [28] $$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Solution in $\Omega_I$: (a)... (d) most representative modes; (e) reduced coordinates evolution; (f) solution at $t = 1$</figcaption></figure>

where

$$ \mathcal{L}(\Psi) = \frac{\partial}{\partial \mathbf{x}} \left( \left[ \nabla \mathbf{v} \cdot \mathbf{x} - \frac{1}{2} f(\mathbf{x}) \cdot \mathbf{x} \right] \cdot \Psi \right) - \frac{1}{2} \frac{\partial^2 \Psi}{\partial \mathbf{x}^2} \quad [29] $$
$$ \mathcal{K}(\Psi) = \partial\Psi/\partial t \qquad [30] $$

and $\nabla v$ is a given velocity gradient of the fluid assumed homogeneous in space.

<figure><img src="image_5.png" /><figcaption>Figure 4. Fully finite element solution in $\Omega_2$</figcaption></figure>

Equation [28] is integrated by assuming that the p.d.f. vanishes on the domain boundary, i.e.

$$ \Psi(\mathbf{x} \in \partial B(0, \sqrt{b}), t) = 0 \qquad [31] $$

as well as a quiescent initial distribution that corresponds to the equilibrium distribution in absence of flow

$$ \Psi(\mathbf{x}, t=0) = \frac{\left(1 - \frac{\mathbf{x}^2}{b}\right)^{b/2}}{\left(2\pi b/2 + b\right)} \qquad [32] $$
In what follows we are solving the just defined model by assuming a given flow whose kinematics is described from its gradient of velocity

$$ \nabla \mathbf{v} = \begin{pmatrix} 0 & 6 \\ 0 & 0 \end{pmatrix} \qquad [33] $$

being the parameter $b=10$. For the calculation of the stress related to the presence of molecules we are computing the virial stress (also known within the rheology community as the Kramer's rule) according to:

$$ \tau(t) = \int_{\Omega} (\mathbf{x} \otimes f(\mathbf{x}) \cdot \mathbf{x}) \Psi(\mathbf{x}, t) d\mathbf{x} \qquad [34] $$

Figure 5 depicts the domain decomposition. In the larger region a reduced modelling is performed, whereas in both small regions a fully finite element discretization is applied. We perform this decomposition because the solution behaviour is known in this case. However, when nothing on the solution is known an error indicator should be integrated in order to define the regions needing for a fully finite element description. The definition of appropriate error indicators constitutes a work in progress.

The computed solution making use of the reduced approximation basis can be described from a finite sums decomposition involving 14 terms:

$$ \Psi(\mathbf{x} \in \Omega_1, t) = \sum_{i=1}^{i=14} \alpha_i F_i(\mathbf{x}) G_i(t) \qquad [35] $$

where the space and time functions were normalized (by dividing by its norm). The four most representative space modes are depicted in Figure 6(a)-(d), where all the functions $G_i(t)$ are also represented. Finally the solution at time $t=5$ is reconstructed and shown in Figure 6(f). Figure 7 depicts the finite element solution in the remaining regions at different times. As in the example described in the previous section, by merging the solutions in both regions at different times, the solution in the whole domain is easily defined. Figure 8 depicts the steady solution in the whole domain and Figure 9 shows the evolution of the extra-stress due to the presence of molecules computed from [34].

In the model considered in this section the mesh contains 951 nodes. Obviously this number corresponds with the number of degrees of freedom if one is using a standard finite element model. However, 777 nodes are located in the region $\Omega_1$ that where reduced to 14 degrees of freedom involved in the reduced approximation. The remaining degrees of freedom (951-777=174) consist of the nodal values at nodes located in region $\Omega_2$. Thus, the size of the final discrete systems scales with the number of nodes located in region $\Omega_2$ because the ones related to the reduced approximation (of the order of 10) are negligible with respect to the former ones.
<figure><img src="image_6.png" /><figcaption>Figure 5. Domain decomposition for the solution of the kinetic theory model</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Solution in $\Omega_I$: (a)… (d) most representative modes; (e) reduced coordinates evolution; (f) solution at $t = 5$</figcaption></figure>
## 5. Conclusions

In this paper we presented a first attempt for coupling reduced and finite element approximation bases. The reduced approximation basis allows significant reduction of the number of degrees of freedom involved in the resulting discrete model in the regions in which the solution evolves smoothly. On the other hand, in the regions where the solution exhibits high gradients or even discontinuities a standard finite element approximation seems better. Coupling of both descriptions is straightforward. Obviously no reduction is expected in the finite element description regions, however, in the remaining one, a reduced number of functions, and then of degrees of freedom. At present the user must specify the domain decomposition, but one can imagine that this task could be easily carried out automatically. The ongoing work consists of coupling separated and FE representations.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Solution in $\Omega_2$: finite element representation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Steady state probability distribution function</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** *Evolution of the extra - shear stress during time*</figcaption></figure>

**6. References**

Ammar A., Ryckelynck R., 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 representation”, *Journal of Non-Newtonian Fluid Mechanics*, 144, 2007, p. 98-121.

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

Holmes P.J., Lumleyc J.L., Berkoozld G., Mattinglya J.C., Wittenberg R.W., “Low-dimensional models of coherent structures in turbulence”, *Physics Reports*, 287, 1997.

Karhunen K., “Uber lineare methoden in der wahrscheinlichkeitsrechnung”, *Ann. Acad. Sci. Fennicae, ser. Al. Math. Phys.*, 37, 1946.

Krysl P., Lall S., Marsden J.E., “Dimensional model reduction in non-linear finite element dynamics of solids and structures”, *International Journal Numerical Methods in Engineering*, 51, 2001, p. 479-504.

Ladeveze P., *Nonlinear computational structural mechanics*, Springer, NY, 1999.
Loève M.M., "Probability theory", *The University Series in Higher Mathematics*, 3rd Ed. Van Nostrand, Princeton, NJ, 1963.

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

A. Nouy, "Recent developments in spectral stochastic methods for the solution of stochastic partial differential equations", *Archives of Computational Methods in Engineering*, vol. 16, n° 3, 2009, p. 251-285.

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

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, 13/1, 2006, p. 91-128.

Sirovich L., "Turbulence and the dynamics of coherent structures part I: Coherent structures", *Quaterly of applied mathematics*, XLV, 1987, p. 561-571.