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

HAL
open science

# On Model Reduction and its Potential Applications in Chemical Engineering

Francisco Chinesta, Amine Ammar, Elías Cueto

► To cite this version:

Francisco Chinesta, Amine Ammar, Elías Cueto. On Model Reduction and its Potential Applications in Chemical Engineering. World Journal of Chemical Engineering, 2007, 1 (1), pp.3-12. hal-01516394v2

HAL Id: hal-01516394

https://hal.science/hal-01516394v2

Submitted on 14 Nov 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 Model Reduction and its Potential Applications in Chemical Engineering

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

¹ LMSP UMR 8106 CNRS-ENSAM-ESEM. 151 Boulevard de l'Hôpital, F-75013 Paris, France

* Corresponding Author: Email: francisco.chinesta@paris.ensam.fr

² Laboratoire de Rhéologie, UMR 5520 CNRS-UJF-INPG. 1301 rue de la piscine, BP 53 Domaine universitaire. 38041 Grenoble Cedex 9, France.

³ Group of Structural Mechanics and Material Modelling. Aragón Institute of Engineering Research (I3A). University of Zaragoza. Maria de Luna, 3. E-50018 Zaragoza, Spain.

Models encountered in chemical engineering usually involve many degrees of freedom related to its multidimensional and multiphysical character. Model reduction is an appealing alternative for treating efficiently (fast and accurately) such models. In this work we will illustrate a few capabilities related to some reduction strategies in order to reduce the complexity of different models. We will analyze the applicability, limitations and enrichment of standard reduction techniques based on the use of the Karhunen-Loève decomposition (also known as proper orthogonal decomposition).

**Keywords:** Model reduction, Karhunen-Loève decomposition, Krylov's subspaces.

## Introduction

Most engineering systems can be modeled by a continuous model, usually expressed as a system of linear and non-linear, partially coupled, 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, meaning the knowledge of the different fields characterizing the physical system at any point and for any 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 only sought at some points and times from which it could be interpolated to any other point and time. Techniques that allow this kind of representation are known as discretization techniques. Today there exist 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 the resolution of complex systems involving millions of unknowns related to the discrete model, possible today. 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 general computation in real time, it is clear that numerous problems must be solved and for this reason the question concerning 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 whom this question has been addressed. For specialists in computational science an answer to this question required the improvement of computational resources, high performance computing, and the use of parallel computing platforms. For some specialists in numerical analysis the challenge lies in the fast resolution of linear systems via the use of preconditioners or multigrid techniques among many others. For others the idea is adapted to the cloud of nodes (points where the solution is computed) in order to avoid an excessive number of unknowns. However, today all these approaches alleviate the computation efforts but the fast and accurate computation remains a real challenge.

This paper presents a different approach, allowing fast and accurate computations based on model reduction. The idea is very simple. Consider a domain where a certain model as well as the associated cloud of nodes is defined. These are able to represent the solution through interpolation. In general the number of unknown scales along with the number of nodes and the solution are smoothly evolving over time. 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 throughout 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). The sizes of the resulting linear problems is very low and CPU time savings in the order of $10^5$ can be attained.

The extraction of this relevant information is a well known topic that is based on the application of the proper orthogonal decomposition, also known as the Karhunen-Loève decomposition (Karhunen, 1946; Loève, 1963) which 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 Park's works, e.g. Park and Cho, 1996). Reduced models usually perform the simulation of some similar problems during a short time interval. The Karhunen-Loève decomposition can be performed as a solution. This will allow for the extraction of the most relevant functions which describe the solution evolution. It is assumed that the solution of a "similar" problem can be expressed using this reduced approximation basis. This reduces the size of the discrete problem and leads to significant CPU time savings. In general, the question concerning 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 then applied at a later stage in different domains (Ryckelynck et al.; 2006, Ammar et al. 2006).

In this paper we illustrate the main ideas of this reduction strategy in a numerical analysis for non specialists. This is done to show its potentiality in the many engineering domains, particularly in the context of chemical engineering.

## Reduced modeling

### The Karhunen-Loève decomposition

We assume that the evolution of a certain field, depending on the physical space $x$ and on the time $t$, $u(x,t)$ is known. In practical applications this field is expressed in a discrete form, meaning that it is known at the nodes of a spatial mesh and for some times $u(x_i, t^p) = u_i^p$. Through the means of writing we can also introduce a spatial Interpolation: $u^p(x) = u(x, t = p\Delta t)$; $\forall p \in [1, \dots, P]$. The main idea behind the Karhunen-Loève (KL) decomposition, is the method of obtaining the most typical or characteristic structure $\varphi(x)$ among these: $u^p(x), \forall p$. This is equivalent to obtaining functions $\varphi(x)$ maximizing $\alpha$

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

which leads to:
$$ \sum_{p=1}^{p=N} \left[ \left[ \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) u^p(\mathbf{x}_i) \right] \left[ \sum_{j=1}^{j=N} \phi(\mathbf{x}_j) u^p(\mathbf{x}_j) \right] \right] = \alpha \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i) \quad (2) $$

∀$\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=N} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \phi(\mathbf{x}_j) \right\} \tilde{\phi}(\mathbf{x}_i) \right] = \alpha \sum_{i=1}^{i=N} \tilde{\phi}(\mathbf{x}_i) \phi(\mathbf{x}_i) \quad (3) $$

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

$$ \tilde{\phi}^T \mathbf{k} \phi = \alpha \tilde{\phi}^T \phi ; \quad \forall \tilde{\phi} \quad \Rightarrow \quad \mathbf{k} \phi = \alpha \phi \quad (4) $$

The two points' correlation matrix is given by

$$ \mathbf{k}_{ij} = \sum_{p=1}^{p=N} u^p(\mathbf{x}_i) u^p(\mathbf{x}_j) \Leftrightarrow \mathbf{k} = \sum_{p=1}^{p=N} \mathbf{u}^p(\mathbf{u}^p)^T \quad (5) $$

which is a 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) $$

It is easy to verify that the matrix **k** in Eq. (5) results in

$$ \mathbf{k} = \mathbf{Q} \mathbf{Q}^T \quad (7) $$

## Reduced modeling

If it is known that the evolution of a certain field

$$ 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) $$

is coming from some direct simulations or from experimental measures, then matrices **Q** and **k** can be computed and the eigenvalue problem given by Eq. (4) solved. The solution of Eq. (4) results in **N** couples eigenvalue-eigenvector. However, in a large number of models involving regular time evolutions of the solution, the magnitude of the eigenvalues decrease 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 > \cdots > \alpha_N$. The $n$ eigenvalues belonging to the interval $\alpha_1 > \cdots > \alpha_n$ with $\alpha_n > \alpha_1 10^{-8}$ and $\alpha_{n+1} < \alpha_1 10^{-8}$ are selected, because their associated eigenvectors are enough to accurately represent the entire solution evolution. In a large variety of models $n \square N$ moreover $n$ only depends on the regularity of the solution evolution, and not on the dimension of the physical space (1D, 2D or 3D) or on the size of the model **N**.

The reduced approximation basis consists of the $n$ eigenvectors $\phi_1, \cdots, \phi_n$, allowing one to define the basis transformation matrix **B**:

$$ \mathbf{B} = (\phi_1, \phi_2, \cdots, \phi_n) \quad (9) $$

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

$$ \mathbf{u} = \sum_{i=1}^{n} \phi_i \xi_i(t) = \mathbf{B} \xi(t) \quad (10) $$

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} \quad (11) $$

where $\mathbf{f}^{p-1}$ accounts for the solution at the previous time step. When we introduce Eq. (10) the result is:

$$ \mathbf{A} \mathbf{u}^p = \mathbf{f}^{p-1} \Rightarrow \mathbf{A} \mathbf{B} \xi^p = \mathbf{f}^{p-1} \quad (12) $$

and when both terms are multiplied by $\mathbf{B}^T$ it results in:

$$ \mathbf{B}^T \mathbf{A} \mathbf{B} \xi^p = \mathbf{B}^T \mathbf{f}^{p-1} \quad (13) $$

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

**Remark 1.** Equation (13) can also be derived by introducing the approximation (10) into the partial differential equation Galerkin form.
## Enriching the reduced approximation basis

The strategy that was just described, allows one to compute large size models very quickly. One could for example solve the full model using a standard discretization technique (finite differences, finite elements ...) for a small time interval. Matrix $Q$ and $k$ can then be defined, computing the reduced approximation basis transformation $B$, which leads to the reduced solution procedure illustrated by Eq. (13).

However, it is not guaranteed that the reduced basis that has been built from the solution, evolving over a short time interval, remains accurate for describing the solution during the entire simulation interval. During the simulation, material properties and 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.

If one would compute reduced model solutions and maintain one's confidence regarding the related solution, an accuracy check must be performed and an enrichment strategy must also be defined in order to adapt the reduced approximation basis. This makes it possible to capture the new events (present in the solution evolutions) which cannot be described accurately from the original reduced approximation basis.

For this purpose, Ryckelynck proposed (Ryckelynck, 2004) that one 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 that has previously been solved, or the ones coming from a full simulation for 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:

$$R = Au^S - f^{S-1} = AB\xi^S - f^{S-1} \quad (14)$$

If the norm of the residual is small enough, $||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, $||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: $B \leftarrow (B, R, AR, A^2R)$.

One could expect that the enrichment process is continuously increasing the size of the reduced approximation basis but in fact, after reaching the convergence, the Karhunen-Loève decomposition is performed on the whole past time interval in order to extract the significant information and to define an orthogonal reduced approximation basis.

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

## Numerical examples

### An introductive example: analyzing the reduced basis capabilities

In this section we are considering a simple 1D model related to the heat transfer equation (we omit the units, all of them being expressed in metric system):

with $k = 0.01$, $t \in ]0,30]$ and $x \in ]0,1[$. The initial condition reads $T(x,t=0)=1$ and the boundary conditions are given by $-k \frac{\partial T}{\partial x}|_{x=0,t} = q(t)$ and $-k \frac{\partial T}{\partial x}|_{x=1,t} = 0$. The boundary source $q(t)$ is prescribed to different values during the simulations that follow.

Equation (15) is discretized by using the implicit finite element method on a mesh that consists of 100 nodes, where a linear approximation is defined in each of the resulting 99 elements. The time step was set to $\Delta t = 0.1$. The resulting discrete system can be written as:

$$KT^p = MT^{p-1} + q^p \quad (16)$$
where vector $q^p$ accounts for the boundary heat flux source at each time step $p$.

Firstly, we are considering the solution of the thermal model as just described. We relate it to the following boundary heat source:

$$ q(t) = \begin{cases} 1 & 0 < t < 10 \\ 0 & t \ge 10 \end{cases} \qquad (17) $$

The computed solution is depicted in figure 1 where the temperature profiles at times $t_n = n$ ($n=1,2,\dots,30$) are represented. The evolution of the temperature profiles within the first 10 seconds (heating stage) is depicted in red. In the remaining time interval no more heating sources exist. The heat is then moving through a conduction mechanism from the more heated zones towards the coldest ones. The profiles within this time interval are represented in blue.

<figure><img src="image_2.png" /><figcaption>Figure 1. Temperature profiles related to the thermal model with the source term modeled by Eq. (17).</figcaption></figure>

Now, from these 30 profiles we could define matrices $Q$ and $k$, which lead to the eigenvalue problem, allowing us to extract the significant eigenvectors according to Eq. (4). The resulting eigenvalues were: $\alpha_1 = 1790$, $\alpha_2 = 1.1$, $\alpha_3 = 0.1$ and $\alpha_j < \alpha_1 \times 10^{-8} (4 \le j \le 100)$.

This result implies that the whole solution evolution could be accurately represented as a linear combination of the 3 eigenvectors related to the first 3 highest eigenvalues. In order to impose the initial condition with ease, it is sometimes suitable to add to these eigenvectors the initial condition (even if the resulting approximation basis is no

more orthogonal). Figure 2 depicts the resulting approximation functions consisting of the 3 eigenfunctions related to the 3 highest eigenvalues along with the initial condition and all of them normalized $\phi_j = \phi_j / \|\phi_j\|$. These functions allow matrix $B$ to be defined as well as the reduced model derived from Eq. (16):

$$ B^T K B \xi^p = B^T (M B \xi^{p-1} + q^p) \qquad (18) $$

which only involves 4 degrees of freedom. Even in the case of non linear models and an implicit discretization, only the inversion of a matrix of size 4 is required at each time step.

If we assume that the initial condition has been placed in the first column of $B$, then the initial condition in the reduced basis writes: $(\xi^0)^T = (1 \ 0 \ 0 \ 0)$. This condition, Eq. (18) can be applied for computing the whole time evolution. According to the basis transformation relationship it is obvious that the global solution can be obtained from the reduced one: $T^p = B\xi^p$. Figure 3 compares some temperature profiles obtained through the use of the global model (Eq. (16)) as well as what was depicted in figure 1 and also the ones obtained from the reduced model (Eq.(18)) (represented by the stars symbol) from which an excellent accuracy can be noticed. This accuracy is not surprising because as indicated before, the 4 approximation functions used were the ones related to the highest eigenvalues which consequently represents the optimal reduced approximation basis.

<figure><img src="image_3.png" /><figcaption>Figure 2. Functions integrating the reduced approximation basis.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 3. Global (continuous line) versus reduced (stars) model solutions.</figcaption></figure>

In order to conclude about the applicability of this reduced approximation basis for simulating problems different to the one that served to compute it, we are considering the thermal model defined in the same domain with the same initial condition but with a slightly different boundary heat source term:

$$ q(t) = \begin{cases} \frac{t}{20} & 0 < t < 20 \\ \frac{t-30}{5} & t \ge 20 \end{cases} \qquad (19) $$

Figure 4 compares the computed reference solution (continuous line) using the model represented by Eq. (16) and the one obtained through the use of the reduced model (18) –stars profiles-. The reduced approximation basis consists of the 4 functions represented in Fig. 2 and that which was associated with the thermal model related to the boundary condition given by Eq. (17). We can notice the excellent accuracy, in some unexpected manner because of the non evident compatibility between the problem, solutions related to Eqs. (17) and (19), and then the unexpected ability of the approximations' functions extracted from the solution of the thermal model defined by Eq. (17) for describing the solution of the thermal model related to Eq. (19).

<figure><img src="image_5.png" /><figcaption>Figure 4. Global (continuous line) versus reduced (stars) model solutions related to the thermal model associated to the thermal source given by Eq. (19).</figcaption></figure>

From this result we start to realize the potentiality of model reduction, but still two questions remain: (i) how do one quantify the quality of a reduced solution without the necessity of computing the global solution? And (ii) in case of noticing a lack of accuracy, how do one enrich the reduced approximation basis in order to improve the solution accuracy?

To address these questions we are considering the technique originally proposed by Ryckelynck (Ryckelynck, 2004) that consists of computing the solution residual defined at a certain time step by step.

$$ \mathbf{R} = \mathbf{K}\mathbf{T}^p - \mathbf{M}\mathbf{T}^{p-1} - \mathbf{q}^p \qquad (20) $$

This residual can be used to quantify the accuracy of the reduced solution and allows one to address the first question. Concerning the second one, we are assuming that the residual resulting from the application of Eq. (20) is greater than the threshold value, i.e. $||\mathbf{R}|| \ge \varepsilon$. A natural election consists of enriching the reduced approximation basis by adding this residual (that is orthogonal in a Galerkin sense to the approximation functions) and some of the Krylov's subspaces related to this residual (according with the procedure previously described) to matrix B.

The residual norm for the thermal model at $t = 30$ was $||\mathbf{R}|| = 2.6 \times 10^{-5}$. This thermal model is related to both the boundary source given by Eq. (19) and the reduced approximation basis depicted in Fig. 2.
This justifies the excellent agreement between the global and reduced solutions noticed in Fig. 4.

Figure 5 depicts the normalized residual $(\mathbf{R} \leftarrow \mathbf{R}/\|\mathbf{R}\|)$. Even if the residual is very low, it can be noticed that the largest deviations are concentrated around the left boundary, where the boundary thermal source applies.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Residual at $t = 30$ related to the thermal solution associated with Eq. (19) and computed using the reduced approximation basis depicted in Fig. 2.</figcaption></figure>

Now, we proceed to enrich the reduced approximation basis by introducing this residual (the first Krylov's subspace) into matrix $\mathbf{B}$ according to: $\mathbf{B} \leftarrow (\mathbf{B}, \mathbf{R})$. If the thermal model related to Eq. (19) is resolved again through the use of the recently updated (enriched) reduced approximation basis, then the norm of the residual decreases. This justifies the introduction of Krylov's subspaces which are related to the residual for the improvement of the reduced solution accuracy. A more in depth analysis on the enrichment procedure is performed in the next section.

### Analyzing the enrichment procedure

We consider in this section the initial domain $\Omega_0 = [-L_0, L_0] \times [-H_0, H_0] = [-\pi, \pi] \times [-\pi, \pi]$ which evolves over time according to the following incompressible velocity field:

$$v(x, y) = \begin{pmatrix} u(x, y) \\ v(x, y) \end{pmatrix} = \begin{pmatrix} \cdot \gamma x \\ \cdot (-\gamma y) \end{pmatrix} \qquad (21)$$

being its length and width at time $t$, $2L_t = 2(L_0 e^{\gamma t})$ and $2H_t = 2(H_0 e^{-\gamma t})$ respectively. Thus, it results in $\Omega_t = [-L_t, L_t] \times [-H_t, H_t]$.

We consider the following linear advection-diffusion problem defined in this domain (the non-linear and coupled case was analyzed in Ryckelynck et al. (2006):

$$\frac{dT}{dt} = k \Delta T \qquad (22)$$

Now, we can proceed to discretize the problem defined by Eq. (22) using a splitting operator technique as well as an implicit Galerkin finite element technique for solving the resulting diffusion problems in the updated geometry. Thus, Eq. (22) is solved in two steps:

In the pure advection explicit step the position of the nodes $\mathbf{x}_i$ is updated according to the nodal velocity, without change in the nodal values:

$$\begin{cases} \mathbf{x}_i^{p+1} = \mathbf{x}_i^p + \mathbf{v}_i^p \Delta t \\ T_i^{p+1/2} (\mathbf{x}_i^{p+1}) = T_i^p (\mathbf{x}_i^p) \end{cases} \qquad (23)$$

In the diffusive step, we solve with an implicit schema the pure diffusion problem in the domain that was just updated according to

$$(\mathbf{M} - \Delta t \mathbf{K}) \mathbf{T}^{p+1} = \mathbf{M} \mathbf{T}^{p+1/2} \qquad (24)$$

We consider the evolution of the temperature field in $\Omega_t$ when a null heat flux is prescribed on the domain boundary $\partial\Omega_t$, i.e. $\nabla T \cdot \mathbf{n}|_{\mathbf{x}\in\partial\Omega_t} = 0$ (being $\mathbf{n}$ the unit outwards vector defined on the domain boundary). The following initial temperature is considered in $\Omega_0$:

$$T(x, y, t=0) = T_0 (1 + \cos(x)) \qquad (25)$$

The numerical simulation is carried out with the choice of parameters grouped in table 1 (all with units in the metric system).

<table><caption>**Table 1:** Simulation parameters considered in the numerical example.</caption><thead><tr><th>Number of nodes: N</th><th>t<sub>max</sub></th><th>L<sub>0</sub></th><th>T<sub>0</sub></th><th>Δt</th><th>• γ</th><th>K</th></tr></thead><tbody><tr><td>101 x 10</td><td>5</td><td>π</td><td>2</td><td>0.1</td><td>0.1</td><td>0.5</td></tr></tbody></table>
It is obvious that the number of nodes assumed in the y-direction is irrelevant due to the unidirectional heat transfer resulting from the initial temperature field (25). This fact facilitates the model validation as well as the solution representation.

Figure 6 depicts each of the temperature profiles in five-time steps. We can notice how the diffusion introduces a smoothing effect at the same time that the domain is growing in the x-direction. We have verified that for a sufficient time-period the temperature becomes uniform, reaching the intended value of the initial temperature given by Eq. (25), i.e. $T(x, y, t \to \infty) \approx T_0$ thus proving the conservative behavior of the discretization schema. Finally, we can also notice the null slope of the temperature at both boundaries, a direct consequence of the boundary condition prescribed.

<figure><img src="image_7.png" /><figcaption>Figure 6. Initial temperature distribution (stars) and temperature profiles (along the x-axis) at times $5\Delta t, 10\Delta t, \dots, 50\Delta t$.</figcaption></figure>

If we denote by $T^{5\Delta t}, \dots, T^{50\Delta t}$ the vectors containing the nodal temperatures at the times $t = 5\Delta t, \dots, 50\Delta t$ respectively, matrix $Q$ reads:

$$Q = \begin{pmatrix} T_1^{5\Delta t} & T_1^{10\Delta t} & \dots & T_1^{50\Delta t} \\ T_2^{5\Delta t} & T_2^{10\Delta t} & \dots & T_2^{50\Delta t} \\ \vdots & \vdots & \ddots & \vdots \\ T_N^{5\Delta t} & T_N^{10\Delta t} & \dots & T_N^{50\Delta t} \end{pmatrix} \quad (26)$$

It is very important to notice that these vectors are associated with different nodal distributions. We solve the eigenvalue problem by extracting the significant functions that need to be taken into consideration in relation to the reduced approximation basis. This result in only two

eigenvalues belonging to the interval defined by the highest eigenvalue and $10^{-8}$ times this value.

Thus, the whole time evolution could accurately be represented as a linear combination of the two associated eigenvectors. We emphasize the fact that both vectors are associated with the nodes but not with a particular nodal distribution. Moreover, in order to represent the initial condition we propose to add the vector form of the initial condition that we denote by $T^0$, to both eigenvectors. As with the example that was covered in the previous section, if this reduced approximation basis consisting of 3 approximation functions is considered, the whole time evolution can be computed with a remarkable accuracy.

Now, we will try to compute a reduced order solution using an adaptive procedure, avoiding any “a priori” knowledge. For this purpose we start with a tentative approximation basis that contains a single function corresponding to the initial condition

$$T(t) = T^0 \quad \zeta(t) = B \xi \quad (27)$$

Introducing the expression of the approximation basis given by Eq. (27) in the equation governing the evolution of the reduced order solution

$$B^T (M - \Delta t K) B \xi^{p+1} = B^T M B \xi^{p+1/2} \quad (28)$$

and taking into account that the verification of the initial condition implies $\xi^0 = 1$, we can compute its evolution, i.e. $\xi^p$, and the associated nodal temperatures $T^p = B\xi^p$. Figure 7 compares the computed solutions, using the reduced basis at times $5\Delta t, 20\Delta t, 35\Delta t$ & $50\Delta t$ with the reference ones.

Now, from the final reduced order solution at time $50\Delta t$ (figure 7, subfigure down-right, curve in red), which we denote by $T(t = 50\Delta t)$, we can compute the residual that allows us to enrich the approximation basis using the Krylov's subspaces.

The residual in our problem (which defines the first Krylov's subspace) is depicted in figure 8. Despite the fact that the residual has been represented in the final configuration (final domain) it will be assumed related to the nodes, and consequently it becomes well defined in each one of the intermediate domains.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Reduced order solutions computed using the approximation basis consisting on the initial condition (red curves) versus de reference ones (blue curves), at times $5\Delta t, 20\Delta t, 35\Delta t, 50\Delta t$ (up-left, up-right, down-left and down-right respectively).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Computed residual at the last time step.</figcaption></figure>

As the residual norm is higher than a sufficient small value, the whole evolution must be recomputed using an enriched basis. Thus, consider the reduced approximation enriched basis obtained by adding the residual just computed: $\mathbf{B} \leftarrow (\mathbf{B}, \mathbf{R})$

Now, the thermal evolution is recomputed in $[0, t_{\max}]$ from $(\xi^0)^T = (1 \ 0)$. The final norm of the residual results in: $\|\mathbf{R}\| = 10^{-10}$, proving that with only two approximation functions the evolution can be accurately represented.

## Conclusion

In this work we describe an algorithm that automatically extracts the characteristic information describing the solution evolution. This technique involves a Karhunen-Loève decomposition leading to an optimal number of approximation functions in tandem with an enrichment strategy which is able to improve the solution accuracy when needed.

Both, the reduction procedure based on the Karhunen-Loève decomposition and the enrichment procedure have been illustrated in some academic examples involving some simple models and geometries but despite the simplicity of the chosen models (for the sake of clarity) the technique works in very complex situations, leading to considerable CPU time savings (some times in the order of hundred thousands).

The main questions concerning this strategy lie in the reduction of mixed formulations coming from the momentum equations in incompressible media —or in general, in the presence of restrictions—. In this case the reduced approximation functions are subjected to the constraint known as the LBB condition. This implies that within the framework of standard finite elements, approximations of a higher order than the pressure one are used for the velocity interpolation. In this manner further
developments will be required in order to enrich the
different approximation subspaces without violating
the LBB stability condition. In some numerical
experiments we have noticed that enriching the
velocity subspace with more functions than the one
associated to the pressure approximation, the LBB
condition is ensured, but there is not a theoretical
proof in this sense. In any case, both the
decomposition and the enrichment should be
compatible with this aforementioned stability
condition and they should be performed without
increasing the computation time.

References

Ammaer, A., Rychelynck, D., Chinesta, F. and
Keunings, R. On the reduction of kinetic theory
models related to finitely extensible dumbbells.
*Journal of Non-Newtonian Fluid Mechanics*,
*134*, 136-147, (2006).

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

Karhunen, K. Über lineare Methoden in der
wahrscheinlichkeitsrechnung. *Ann. Acad. Sci.
Fennicae, ser. Al. Math. Phys.*, 37 (1946).

Krysl, P., Lall, S., and Marsden, J.E. Dimensional
model reduction in non-linear finite element
dynamics of solids and structures. *Int. J.
Numer. Meth. in Engng.*, 51, 479–504 (2001).

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

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

Ryckelynck, D. A Priori Hyperreduction Method: an
Adaptive Approach. *Journal of Computational
Physics*, 202, 346-366 (2005).

Ryckelynck, D., Chinesta, F., Cueto, E. and
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, 91-128 (2006).

Sirovich, L. Turbulence and the dynamics of coherent structures part I: Coherent structures. *Quaterly of applied mathematics*, XLV, 561–571 (1987).