# On the a priori Model Reduction: Overview and Recent Developments

David Ryckelynck, Francisco Chinesta, Elías Cueto, Amine Ammar

► To cite this version:

David Ryckelynck, Francisco Chinesta, Elías Cueto, Amine Ammar. On the a priori Model Reduction: Overview and Recent Developments. Meshfree Methods for Partial Differential Equations II, 43, 2005, Lecture Notes in Computational Science and Engineering, 10.1007/BF02905932. hal-01633395

HAL Id: hal-01633395

https://hal.science/hal-01633395

Submitted on 12 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the *a priori* Model Reduction: Overview and Recent Developments

D. Ryckelynck and F. Chinesta

LMSP UMR 8106 CNRS-ENSAM-ESEM
151 Boulevard de l'Hôpital, F-75013 Paris, France
david.ryckelynck@paris.ensam.fr

E. Cueto

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

A. Ammar

Laboratoire de Rhéologie
1301 rue de la piscine, BP 53 Domaine universitaire
38041 Grenoble Cedex 9, France

## Summary

Karhunen-Loève expansion and snapshot POD are based on principal component analysis of series of data. They provide basis vectors of the subspace spanned by the data. All the data must be taken into account to find the basis vectors. These methods are not convenient for any improvement of the basis vectors when new data are added into the data base. We consider the data as a state evolution and we propose an incremental algorithm to build basis functions for the decomposition of this state evolution. The proposed algorithm is based on the APHR method (A Priori Hyper-Reduction method). This is an adaptive strategy to build reduced order model when the state evolution is implicitly defined by non-linear governing equations. In case of known state evolutions the APHR method is an incremental Karhunen-Loève decomposition. This approach is very convenient to expand the subspace spanned by the basis functions. In the first part of the present paper the main concepts related to the “a priori” model reduction technique are revisited, as a previous task to its application in the cases considered in the next sections.

Some engineering problems are defined in domains that evolve in time. When this evolution is large the present and the reference configurations differ significantly. Thus, when the problem is formulated in the total Lagrangian framework frequent remeshing is required to avoid too large distortions of the finite element mesh. Other possibility for describing these models lies in the use of an updated formulation in which the mesh is conformed to each intermediate configuration. When the finite element method is used, then frequent remeshing must be carried out to perform an optimal meshing at each intermediate configuration. However, when the natural element method, a novel meshless technique, is considered, whose accuracy does not depend significantly on the relative position of the nodes, then large simulations can be performed without any remeshing stage, being the nodal position at each intermediate configuration defined by the transport of the nodes by the material velocity or the advection terms. Thus, we analyze the extension of the “a priori” model reduction, based on the use in tandem of the Karhunen-Loève decomposition (that extracts significant information) and an approximation basis enrichment based on the use of the Krylov’s subspaces, previously proposed in the framework of fixed mesh simulation, to problems defined in domains evolving in time.

Finally, for illustrating the technique capabilities, the “a priori” model reduction will be applied for solving the kinetic theory model which governs the orientation of the fibers immersed in a Newtonian flow.

# 1 THE APHR METHOD: AN INCREMENTAL KARHUNEN-LOÈVE EXPANSION

## 1.1 Introduction

This work focuses in the development of a method for the reduction of the number of state variables of non-linear thermomechanical time dependent problems. This is the A Priori Hyper-Reduction method (APHR) [14]. The main aim of this section is to explain how this
method can be used as an incremental and adaptive version of the Karhunen-Loève expansion. The Karhunen-Loève expansion, [8] [10] (KL expansion) and the Proper Orthogonal Decomposition, [11, 7, 6] (POD) are two names for the same method. This method provides basis functions for a simplified representation of state evolutions, or set of data. These basis functions are eigenfunctions of the averaged autocorrelation of known state evolutions. It can be shown that the KL expansion yields an optimal set of orthonormal basis functions in the sense that the fewest number of functions of all possible bases are required for a given level of accuracy in reconstructing the original set of data [2]. Those basis functions are usually called the empirical eigenvectors. This method can be applied in the framework of signal processing, experimental analysis or computational method in order to build reduced order models. The empirical eigenvectors corresponding to the highest eigenvalues are the functions that describe the main significant events involved in the state evolution. The size of the eigenproblem is equal to the number of state variables. So, in case of large number of state variables the eigenproblem becomes high time consuming. To overcome this difficulty Sirovich proposed the snapshot POD [15]. Approximate empirical eigenvectors are obtained by a linear combination of few selected states at different time instants: the snapshots. The coefficients of this linear combination are obtained thanks to an eigenproblem as large as the number of states selected to perform the decomposition. But, in case of transient state evolutions the snapshot POD can not reproduce the events having a characteristic time lower than the time step between snapshots. So it is difficult to choose this time step properly. The basis functions used for the decomposition of the state evolution depend on the evolutions considered to perform the KL expansion or the snapshot POD. In the scientific community of signal processing this approach is considered as an adaptive one while, for example, the basis functions of the Discret Consin Transform do not depend on what the state evolution could be. In this article we do not address this meaning of the adaptive term. The main drawback of the KL expansion and snapshot POD is that no adaptation of the empirical eigenvectors is possible when a new state evolution has to be considered in order to improve the basis functions or in order to increase the number of basis functions [9]. A full snapshot POD or KL expansion must be done.

The purpose of the adaptive method that we propose is to simplify the computation of the new basis functions by taking into account a first set of basis functions and the new state evolution. To do so we propose an incremental approach: the APHR method. This method was developed for the decomposition of a state evolution implicitly defined by a set of non-linear time dependent equations. This implicit case leads to a model reduction method [14]. When the state evolution is known, the APHR method leads to an incremental Karhunen-Loève transform. A simple example is introduced to illustrate the incremental Karhunen-Loève transform.

## 1.2 The APHR Method

Let's consider the evolution of a state field $s_{ref}(M, t)$ defined at any time $t$ of the time interval $[0, T]$ and defined for any point $M$ in a domain $\Omega$. The purpose of the KL expansion of $s$ is to find a set of basis functions $(\Psi_k(M))_{k=1...r}$ defined in $\Omega$ and a set of reduced state variables $(v_k(t))_{k=1...r}$ defined in $[0, T]$ such that:

1. Preliminary stage: $\Psi_k$ is defined by the stationarity of $\lambda(\Psi)$

$$ \lambda(\Psi) = \frac{\int_{0}^{T} \left( \int_{\Omega} s_{ref}(M, t) \Psi(M) d\Omega \right)^2 dt}{\int_{\Omega} \Psi^2(M) d\Omega} \quad (1) $$

$$ \delta\lambda(\Psi_k) = 0 \quad (2) $$
$$ \lambda(\Psi_k) \neq 0 \qquad (3) $$

$$ \int_{\Omega} \Psi_k^2(M) \, d\Omega = 1 \qquad (4) $$

2. Projection stage: $(v_k(t))_{k=1...r}$ minimize $\eta(v_1, ..., v_r)$

$$ \eta(v_1, ..., v_r) = \int_0^T \int_{\Omega} \left( s_{ref}(M, t) - \sum_{k=1}^{k=r} \Psi_k(M) v_k(t) \right)^2 d\Omega \, dt \qquad (5) $$

Obviously, this method can be applied on numerical functions, if we choose numerical integration scheme and shape functions $(N_i)_{i=1,...,N}$ ($N$ being the number of nodes) such that:

$$ s_{ref}(M, t) = \sum_{i=1}^{i=N} N_i(M) q_i(t) \qquad (6) $$

The adaptive method that we propose allows to avoid the preliminary construction of basis functions before the incremental decomposition of the state evolution. Obviously, an initial set of known basis functions can be used, but it is not necessary. During the decomposition of the state evolution, approximate basis functions $(\phi_k^{(m)})_{k=1...r^{(m)}}$ are improved. The superscript $m$ indicates that the approximation basis has been updated $m$ times until now. We will prove later that the basis functions $\phi_k^{(m)}$ are approximations of the empirical eigenfunction $\Psi_k$. An approximate state evolution $s_{\phi}^{(m)}$ is given by:

$$ s_{\phi}^{(m)}(M, t) = \sum_{k=1}^{k=r^{(m)}} \phi_k^{(m)}(M) a_k^{(m)}(t) \qquad (7) $$

The adaptive procedure includes extensions of the subspace spanned by the basis functions and selections of the most significant basis functions in order to represent the state evolution.

The basis functions being defined in $\Omega$, any linear combination of these basis functions will be defined everywhere in $\Omega$. We expect that the linear combination enables to estimate the state anywhere in the domain $\Omega$. So the reduced state variables can be computed by considering only a part of the whole domain: a reduced integration domain $\Omega_{hyp}$. It is therefore possible to reduce the number of integration points used to compute the reduced state variables. This strategy establishes a model hyperreduction.

Let's consider a weak form of the equations defining the state evolution. The state evolution we want to forecast is $s$. At any time instant, the field $s$ belongs to a state space $S$. We consider a test function $s^*$ in $S$ to obtain the following governing equation:

$$ s(M, 0) = s_{ini}(M) \qquad (8) $$

$$ R(s^*, s, \frac{\partial s}{\partial t}, t) = 0 \quad \forall s^* \in S \qquad (9) $$

$R$ is the residual of the governing equations. Let $s_{ref}$ be the exact solution of these governing equations. The norm of the residual $R(s^*, s_{\phi}^{(m)}, \frac{\partial s_{\phi}^{(m)}}{\partial t}, t)$ provides an error estimator to
check the quality of $s_{\phi}^{(m)}$. A function decomposition can be considered appropriate when this estimator is lower than a given bound $\epsilon_R$ (small enough) for all time $t$ in $[0, T]$. When $s_{ref}$ is known $R$ is defined by:

$$R(s^*, s_{\phi}^{(m)}, t) = \int_{\Omega} s^*(M) (s_{\phi}^{(m)}(M, t) - s_{ref}(M, t)) d\Omega \quad (10)$$

Different adaptive strategies can be developed with the APHR method. The common items of such strategies are:

1. A set of basis functions always defined in the whole time interval;

2. An extension of the subspace spanned by the basis functions thanks to residuals of the governing equations;

3. A selection of basis functions based on a Karhunen-Loève expansion of the evolution of reduced state variables;

4. A formulation allowing compute the reduced state variables with a reduced integration domain $\Omega_{hyp}$ smaller than the whole domain $\Omega$.

The adaptive strategy must precise: (i) when the reduced state variables should be computed; (ii) when the quality of the approximate state should be checked and (iii) when the basis functions should be adapted. The strategy we propose is an incremental approach that proceeds as we describe in the next paragraphs. Let's assume that the decomposition was done properly over the time interval $[0, t_{\alpha}]$ thanks to reduced state variables $a_k^{(m)}$ associated to the basis functions $\phi_k^{(m)}$. The first step of the algorithm is to find the next time instant $t_{\beta}$ ($t_{\beta} \ge t_{\alpha}$) such that the residual computed at $t_{\beta}$ does not satisfy the quality criteria. When $t_{\beta}$ is known, an adaptation of the basis functions is performed until the quality criteria is satisfy at $t_{\beta}$. When the decomposition at $t_{\beta}$ is good enough, the incremental computation-adaptation can be continued. If the initial basis functions $\phi_k^{(0)}$ does not allow a convenient decomposition of $s_{ini}$ the first $t_{\beta}$ is 0. Obviously, if no basis function is known, we can choose $\phi_1^{(0)}$ collinear to $s_{ini}$.

For the case where $s_{ref}$ is known, the adaptation is done in two steps to satisfy the quality criteria:

1. The first step is the KL expansion of the reduced state variables over the time interval $[0, t_{\beta}]$. We denote $\underline{V}_j^{(m)}$ the eigenvectors of the matrix $\underline{C}^{(m)}$ of the averaged autocorrelations such that:

$$C_{kp}^{(m)} = \int_{0}^{t_{\beta}} a_{k}^{(m)}(t) a_{p}^{(m)}(t) dt \quad k, p \in \{1, \dots, r^{(m)}\} \quad (11)$$

Let $\mu_j^{(m)}$ be the eigenvalue related to the eigenvector $\underline{V}_j^{(m)}$ such that:

$$\mu_1^{(m)} \ge \mu_2^{(m)} \ge \dots \ge \mu_{r^{(m)}}^{(m)} \quad (12)$$

with

$$\underline{V}_i^{(m)} \underline{V}_j^{(m)} = \sum_{k=1}^{k=r^{(m)}} \underline{V}_k^{(m)} \underline{V}_{kj}^{(m)} = \delta_{ij} \quad (13)$$
where $V_k$ refers the k-component of vector $\underline{V}_i$

The selection of the most significant basis functions is performed by taking the $\tilde{r}^{(m+1)}$ eigenvectors related to the eigenvalues greater than $\mu_{\tilde{r}^{(m+1)}}^{(m)}$, with:

$$ \mu_{\tilde{r}^{(m+1)}}^{(m)} \geq \varepsilon_{KL} \mu_{1}^{(m)} \quad (14) $$

where $\varepsilon_{KL}$ is a small enough parameter (usually we choose $\varepsilon_{KL} = 10^{-8}$). The selected eigenvectors define the basis reduction matrix $\underline{V}^{(m)}$:

$$ \underline{V}^{(m)} = [\underline{V}_{1}^{(m)}, \dots, \underline{V}_{\tilde{r}^{(m+1)}}^{(m)}] \quad (15) $$

The basis reduction matrix provides new basis functions $(\phi_k^{(m+1)})_{k=1...\tilde{r}^{(m+1)}}$ such that:

$$ \phi_k^{(m+1)}(M) = \gamma_k \sum_{p=1}^{p=r^{(m)}} \phi_p^{(m)}(M) V_{pk}^{(m)} \quad k = 1, \dots, \tilde{r}^{(m+1)} \quad (16) $$

where $\gamma_k$ is computed in order to normalize the basis functions, i.e.:

$$ \int_{\Omega} \phi_k^{(m+1)}(M) \phi_k^{(m+1)}(M) d\Omega = 1 \quad (17) $$

2. The second step of the adaptation is the extension of the subspace spanned by the functions $(\phi_k^{(m+1)})_{k=1...\tilde{r}^{(m+1)}}$ by adding the function $\phi_{\tilde{r}^{(m+1)}}^{(m+1)}(r^{(m+1)} = \tilde{r}^{(m+1)} + 1)$ collinear to the residual at $t_\beta$:

$$ \int_{\Omega} s^{*}(M) \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(M) d\Omega = \gamma_{\tilde{r}^{(m+1)}} R(s^{*}, s_{\phi}^{(m)}, t) \quad \forall s^{*} \in S, t = t_{\beta} \quad (18) $$

where $\gamma_{\tilde{r}^{(m+1)}}$ is chosen to verify:

$$ \int_{\Omega} \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(M) \phi_{\tilde{r}^{(m+1)}}^{(m+1)}(M) d\Omega = 1 \quad (19) $$

Because of the basis adaptation, the reduced state variables have to be updated over $[0, t_\beta]$ according to:

$$ a_k^{(m+1)}(t) = \frac{1}{\gamma_k} \sum_{p=1}^{p=r^{(m)}} a_p^{(m)}(t) V_{pk}^{(m)} \quad \forall t < t_\beta, \forall k < r^{(m+1)} \quad (20) $$

and

$$ a_{\tilde{r}^{(m+1)}}^{(m+1)}(t) = 0 \quad \forall t < t_{\beta} \quad (21) $$
Equation (20) can be derived easily from

$$ \int_{\Omega} \left\{ \sum_{k=1}^{k=\tilde{r}^{(m+1)}} \frac{\phi_l^{(m+1)}(M)}{\gamma_l} \phi_k^{(m+1)}(M) a_k^{(m+1)} \right\} d\Omega = \\ = \int_{\Omega} \left\{ \sum_{p=1}^{p=r^{(m)}} \frac{\phi_l^{(m+1)}(M)}{\gamma_l} \phi_p^{(m)}(M) a_p^{(m)} \right\} d\Omega; \quad \forall l \in [1, \dots, \tilde{r}^{(m+1)}] \quad (22) $$

Now, using the matrix form of Eq.(16):

$$ \underline{\phi}^{(m+1)T}(M) = \underline{\phi}^{(m)T}(M) \underline{V} \underline{\gamma} \quad (23) $$

(with $\underline{\gamma}$ the diagonal matrix containing the $\gamma_k$ coefficients), Eq.(22) can be written in the following matrix form:

$$ \int_{\Omega} \left\{ \underline{V}^T \underline{\phi}^{(m)}(M) \underline{\phi}^{(m)T}(M) \underline{V} \underline{\gamma} \underline{a}^{(m+1)} \right\} d\Omega = \\ = \int_{\Omega} \left\{ \underline{V}^T \underline{\phi}^{(m)}(M) \underline{\phi}^{(m)T}(M) \underline{a}^{(m)} \right\} d\Omega \quad (24) $$

and taking into account the orthogonality of the approximation functions (proved in the next theorem), that establishes

$$ \int_{\Omega} \underline{\phi}^{(m)}(M) \underline{\phi}^{(m)T}(M) d\Omega = \underline{I} \quad (25) $$

with $\underline{I}$ the unit tensor, as well as the orthogonality of vectors $\underline{V}_i$, i.e. $\underline{V}^T \underline{V} = \underline{I}$, it results:

$$ \underline{a}^{(m+1)}(M) = \gamma^{-1} \underline{V}^T \underline{a}^{(m)} \quad (26) $$

that is in fact the matrix form of Eq.(20) □.

For $t \ge t_\beta$, the reduced state variables are computed thanks to a modified Galerkin procedure by introducing a reduced integration domain $\Omega_{hyp}$ such that:

$$ s^*(M) = \sum_{k=1}^{k=r^{(m+1)}} \phi_k^{(m+1)}(M) a_k^* \quad \forall M \in \Omega_{hyp} \quad (27) $$

$$ s^*(M) = 0 \quad \forall M \in \Omega - \Omega_{hyp} \quad (28) $$

$$ R(s^*, s_{\phi}^{(m+1)}, t) = 0 \quad \forall s^*, t > t_{\beta} \quad (29) $$

The construction of $\Omega_{hyp}$ depends on each particular problem. In this work we just show that it is possible to perform an hyperreduction. More details on the choice of $\Omega_{hyp}$ can be found in [14].

To get optimal basis functions, a last KL expansion of the reduced state variables can be performed over $[0, T]$ at the end of the incremental computation.
**Theorem 1.1.** In the case of $s_{ref}$ known, if $\Omega_{hyp} = \Omega$, then the basis functions are orthogonal.

*Proof.* Let's assume that $(\phi_k^{(m)})_{k=1...r^{(m)}}$ are the orthogonal basis functions verifying:

$$ \int_{\Omega} \phi_i^{(m)}(M) \phi_j^{(m)}(M) d\Omega = \delta_{ij} \quad (30) $$

Whatever the selection criteria, the KL expansion of the reduced state variables $(a_k^{(m)})_{k=1...r^{(m)}}$ provides orthogonal basis functions $(\phi_k^{(m+1)})_{k=1...\tilde{r}^{(m+1)}}$:

$$ \int_{\Omega} \phi_i^{(m+1)}(M) \phi_j^{(m+1)}(M) d\Omega = \sum_{k=r^{(m)}}^{k=r^{(m+1)}} \sum_{p=r^{(m)}}^{p=r^{(m+1)}} \gamma_{ij} V_{ki}^{(m)} V_{pj}^{(m)} \int_{\Omega} \phi_k^{(m)}(M) \phi_p^{(m)}(M) d\Omega \quad \forall i \le \tilde{r}^{(m+1)} \forall j \le \tilde{r}^{(m+1)} \quad (31) $$

Taking into account Eq.(30):

$$ \int_{\Omega} \phi_i^{(m+1)}(M) \phi_j^{(m+1)}(M) d\Omega = \gamma_i \gamma_j \sum_{k=1}^{k=r^{(m)}} V_{ki}^{(m)} V_{kj}^{(m)} = \gamma_i \gamma_j \delta_{ij} \quad \forall i \le \tilde{r}^{(m+1)} \forall j \le \tilde{r}^{(m+1)} \quad (32) $$

Finally, for $j = r^{(m+1)}$, the basis function is collinear to the residual which is orthogonal to the rest of the basis functions if a classical Galerkin procedure is used to compute the reduced state variables when $\Omega_{hyp} = \Omega$. So the adaptation of the basis functions provides orthogonal functions. □

**Theorem 1.2.** In the case of $s_{ref}$ known, when a full integration is considered ($\Omega_{hyp} = \Omega$), the decomposition obtained at the end of the incremental computation is the KL decomposition of the approximate state evolution $s_{\phi}^{(m)}$. So, it is an approximate KL expansion of $s_{ref}$ according to the precision $\epsilon_R$.

*Proof.* The empirical eigenfunction $\Psi(M)$ of the KL expansion of $s_{\phi}^{(m)}$ maximizes the projection $\lambda(\Psi)$:

$$ \lambda(\Psi) = \frac{\int_{0}^{T} \left( \int_{\Omega} s_{\phi}^{(m)}(M,t) \Psi(M) d\Omega \right)^2 dt}{\int_{\Omega} \Psi^2(M) d\Omega} \quad (33) $$

Let $\Psi_p$ be the $p^{th}$ empirical eigenvector and $\lambda_p$ its related eigenvalue. Since $(\phi_k^{(m)})_{k=1...r^{(m)}}$ is a basis of a subspace of $S$, we can find $(b_k p)_{k=1...r^{(m)}}$ and $\xi_p$ such that:

$$ \Psi_p(M) = \sum_{k=1}^{k=r^{(m)}} \phi_k^{(m)}(M) b_k p + \xi_p(M) \quad (34) $$

$$ \int_{\Omega} \phi_k^{(m)}(M) \xi_p(M) d\Omega = 0 \quad \forall k \in \{1, \dots, r^{(m)}\} \quad (35) $$
Due to the orthogonality of the basis functions $(\phi_k^{(m)})_{k=1...r^{(m)}}$, the eigenproblem defined by the stationarity of $\lambda(\Psi_p)$ results in the following two equations:

$$ \lambda_p b_{ip} = \sum_{j=1}^{j=r^{(m)}} \int_0^T a_i(t) a_j(t) dt b_{jp} \quad (36) $$

and

$$ \lambda_p \int_{\Omega} s^*(M) \xi_p(M) d\Omega = \int_0^T \int_{\Omega} \int_{\Omega} s^*(M) s_{\phi}^{(m)}(M,t) s_{\phi}^{(m)}(M',t) \xi_p(M') d\Omega d\Omega' dt \quad \forall s^* \in S \quad (37) $$

As

$$ \int_{\Omega} s_{\phi}^{(m)}(M', t) \xi_p(M') d\Omega' = \sum_{k=1}^{k=r^{(m)}} a_k(t) \int_{\Omega} \phi_k^{(m)}(M') \xi(M') d\Omega' = 0 \quad \forall t \quad (38) $$

Eq. (37) implies that $\xi_p$ is equal to zero or it is an eigenfunction associated with $\lambda_p = 0$.

Moreover due to the last KL expansion of the reduced state variables:

$$ \int_0^T a_i^{(m)}(t) a_j^{(m)}(t) dt = \mu_i^{(m)} \delta_{ij} \quad (39) $$

So, the basis functions $(\phi_k^{(m)})_{k=1...r^{(m)}}$ are the empirical eigenfunctions of the Karhunen-Loève expansion of $s_{\phi}^{(m)}$, and the eigenvalue related to $\phi_k^{(m)}$ is $\mu_k^{(m)}$. $\square$

## 2 EXAMPLE OF AN INCREMENTAL KARHUNEN-LOÈVE EXPANSION: A FORMAL DESCRIPTION

In this section we present a formal example of incremental Karhunen-Loève expansion of a given state evolution. The incremental decomposition is started with one initial basis function. Formal computations have been chosen in order to avoid any artefact related to both the numerical integration and the numerical discretisation.

Let's consider the function $s_{ref}(x, t) = \frac{x}{100} + \sin(x) t + \frac{x^2}{(0.1+t)}$ (Figure 1) defined in $\Omega = [-100, 100]$, $\forall t \in [0, T]$ with $T = 40$. To show the capability of the hyperreduction, we choose $\Omega_{hyp} = [-1, 5]$ (Figure 1). A dichotomic search of $t_{\beta}$ has been used in this example. Another approach for finding $t_{\beta}$ can be found in [14].

We also choose as a first basis function a constant field $\phi_1^{(o)}(x) = \frac{1}{\sqrt{200}} (r^{(o)} = 1)$. Moreover we have considered the following algorithm parameters: $\varepsilon_{KL} = 10^{-8}$, $\varepsilon_R = 10^{-3} \|s_{ref}\|_{\Omega}$ (with $\|s_{ref}\|_{\Omega}^2 = \int_{\Omega} s_{ref}^2 d\Omega$). The different steps of the incremental Karhunen-Loève decomposition of $s_{ref}$ are:

* The first $t_{\beta}$ is 0. We obtain: $a_1^{(o)}(0) = 990.232$

* First adaptation:

- Karhunen-Loève expansion of $a_1^{(o)}$ in $[0, 0]$: $\phi_1^{(1)} = \phi_1^{(o)}$
<figure><img src="image_2.png" /><figcaption>Figure 1. $s_{ref}(x, t)$ over $\Omega \times [0, T]$ and a zoom over $\Omega_{hyp} \times [0, T]$</figcaption></figure>

- Extension, $r^{(1)} = 2$: $\phi_2^{(1)}(x) = 1.58299 \cdot 10^{-5}(-2.64563 + x)(2.64663 + x)$ (residual at $t = 0$).

* Computation of the reduced state variables considering $\Omega_{hyp}$ instead of $\Omega$, for $t \ge 0$:

$$a_1^{(1)}(t) = 0.28284 + 0.60491 t + \frac{98.995}{0.1 + t}$$

$$a_2^{(1)}(t) = \frac{-3839}{0.1 + t} (-4.0245 + t)(4.0888 + t)$$

* The state evolution is not admissible for $t > 0.1$, because

$$\|R\|_{\Omega} = 0.000996 \quad \|s_{ref}\|_{\Omega} \quad \text{at } t = 0.1$$

Thus, we set $t_{\beta} = 0.1$.

* Second adaptation:

- Karhunen-Loève expansion of $a_k^{(1)}$ in $[0, 0.1]$, and selection of the unique significant eigenvector:

$$\underline{V}^{(1)} = \begin{Bmatrix} 0.0015682 \\ 0.99999 \end{Bmatrix}$$

$$\phi_1^{(2)}(x) = 4.8363 \cdot 10^{-8} + 1.5811 \cdot 10^{-8} x + 1.5811 \cdot 10^{-5} x^2$$

- Extension from the residual at $t = 0.1$, $r^{(2)} = 2$:

$$\phi_2^{(2)}(x) = 1.5818 \cdot 10^{-5} (-1.3187 + x)(2.3217 + x) + 3.17 \cdot 10^{-4} \sin(x)$$

- The reduced state variables are updated in $[0, 0.1[$

$$a_1^{(2)}(t) = -\frac{3843.4 (-4.0246 + t)(4.0888 + t)}{0.1 + 1.1t}$$

$$a_2^{(2)}(t) = 0$$
* Computation of the reduced state variables considering $\Omega_{hyp}$ instead of $\Omega$, for $t \ge 0.1$:

$$ 
\begin{aligned}
a_1^{(2)}(t) &= -\frac{3323.6(-4.3333 + t)(4.3906 + t)}{0.1 + t} \\
a_2^{(2)}(t) &= \frac{3217.1(-0.011866 + t)(0.12253 + t)}{0.1 + t}
\end{aligned}
$$

* The state evolution is not admissible for $t > 4.9$, because

$$ \|R\|_{\Omega} = 0.00102 \|s_{ref}\|_{\Omega} \quad \text{at } t = 4.9 $$

Thus, we set $t_{\beta} = 4.9$.

* Third adaptation:

– Karhunen-Loève expansion of $a_k^{(2)}$ in [0, 4.9] which results in two significant eigenvectors:

$$ \underline{V}^{(2)} = \begin{Bmatrix} 0.999909 & -0.0135255 \\ 0.0135255 & 0.999909 \end{Bmatrix} $$

$$ \phi_1^{(3)}(x) = 1.5812 \ 10^{-5} (-0.18752 + x) (0.2019 + x) + 4.2305 \ 10^{-6} \sin(x) $$

$$ \phi_2^{(3)}(x) = 1.5818 \ 10^{-5} (-1.3253 + x) (2.342 + x) + 3.2133 \ 10^{-4} \sin(x) $$

– Extension from the residual at $t = 4.9$, $r^{(3)} = 3$:

$$ \phi_3^{(3)}(x) = 1.3507 \ 10^{-5} (-40.556 + x) (-3.4517 + x) - 2.6478 \ 10^{-4} \sin(x) $$

– The reduced state variables are updated in [0, 4.9], being $a_3^{(3)} = 0$.

* Computation of the reduced state variables considering $\Omega_{hyp}$ instead of $\Omega$ for $t \ge 4.9$.

* The state evolution is not convenient after $t = 10.5$.

$$ \|R\|_{\Omega} = 0.00102 \|s_{ref}\|_{\Omega} \quad \text{at } t = 10.5 $$

* Fourth adaptation:

– Karhunen-Loève expansion of $a_k^{(3)}$ in [0, 10.5] which results in two significant eigenvectors:

$$ \underline{V}^{(3)} = \begin{Bmatrix} 0.997634 & -0.0687428 \\ -0.068728 & -0.997298 \\ -0.00147 & -0.0259241 \end{Bmatrix} $$

$$ \phi_1^{(4)}(x) = 1.5811 \ 10^{-5} (-0.01327 + x) (0.012948 + x) - 1.8838 \ 10^{-5} \sin(x) $$

$$ \phi_2^{(4)}(x) = -1.2621 \ 10^{-8} - 5.924 \ 10^{-7} x - 1.5811 \ 10^{-5} x^2 - 2.8834 \ 10^{-4} \sin(x) $$

– Extension from the residual at $t = 10.5$, $r^{(4)} = 3$:

$$ \phi_3^{(4)}(x) = 1.5796 \ 10^{-5} (-3.5587 + x) (0.081079 + x) + 6.0791 \ 10^{-5} \sin(x) $$
- The reduced state variables are updated over $[0, 10.5[$, being $a_3^{(4)} = 0$

* Computation of the reduced state variables considering $\Omega_{hyp}$ instead of $\Omega$ for $t \ge 10.5$.

* The state evolution is convenient until $t = 40$

We obtain a very accurate decomposition :

$$ \sqrt{\frac{\int_0^T \int_{\Omega} (s_{\phi}^{(4)} - s_{ref})^2 d\Omega dt}{\int_0^T \int_{\Omega} s_{ref}^2 d\Omega dt}} = 2 \cdot 10^{-5} $$

As for the snapshot POD, different time instants are selected to expand the subspace spanned by the basis functions. But in the proposed approach they are selected according to the value of the residual due to the current basis functions. Moreover, since each extension is obtained with the residual we are sure to improve the subspace spanned by the basis functions.

A snapshot POD of $s_{\phi}^{(4)}$ can be performed when the decomposition is done. We can choose as snapshots the states corresponding to the time instants of the adaptations: $\{0, 0.1, 4.9, 10.5\}$. The eigenvalues provided by the snapshot POD are:

$$ \{5.00195 \cdot 10^{11}, 13543.8, 25.741, 2.38875 \cdot 10^{-5}\} $$

The eigenvalues of the autocovariance matrix of the reduced state variables, for the whole time interval are:

$$ \{4.55512 \cdot 10^{11}, 1.962379 \cdot 10^{10}, 249052\} $$

For both approaches there are only 3 significant basis functions for the same accuracy. The snapshot POD provides eigenvectors that enable to reduce the contribution of the third empirical eigenfunction. It also improves the contribution of the first basis function. If we choose a full integration scheme ($\Omega_{hyp} = \Omega$) during the incremental Karhunen-Loève transform, the Karhunen Loève decomposition of the reduced state variables over $[0, T]$ provides the following eigenvalues:

$$ \{5.00195 \cdot 10^{11}, 13538.9, 3.86629 \cdot 10^{-5}\} $$

This decomposition is better than the one obtained with the snapshot POD at the end of the first incremental decomposition because it reduces the contribution of the third empirical eigenfunction.

# 3 FIRST ATTEMPTS IN COUPLING MODEL REDUCTION AND MESHLESS TECHNIQUES

## 3.1 The Problematic of Evolving Domains and Discretisations Involving a High Number of Degrees of Freedom

Some engineering problems are defined in domains that evolve in time. The consideration of a moving mesh strategy (as used for example in the Lagrangian finite element formulations) allows to describe accurately the domain evolution, although one must be careful in the boundary tracking due to the eventual contacts for example. Moreover, the advection terms can be accurately integrated using the method of characteristics along the nodal trajectories. However, as it is well known in the context of the Lagrangian finite element method, the mesh becomes too distorted in few iterations to guarantee an accurate field
interpolation in the mesh elements. In order to alleviate the remeshing constraint, some
meshless methods have been proposed. However, usual meshless techniques do not define a
nodal interpolation, and in consequence important difficulties are found in the imposition of
the essential boundary conditions. The Natural Element Method —NEM— [16], is a novel
meshless method, which has the property of nodal interpolation, and its accuracy does
not depend on the regularity of nodal distribution, i.e. there is not geometrical restriction
in the relative position of the nodes. Thus, if the NEM is used in the discretisation of
the variational formulation of the problem, then, the nodal position can be updated using
the material velocity or the advection field, at the same time that advection equations are
integrated using the method of characteristics. Even in the case of very irregular nodal
distributions, when the solution can be interpolated by using the approximation functional
basis, no remeshing is required. In any case, the introduction or the elimination of some
nodes is an easy task [12].

Some 2D or 3D engineering problems remain today untreatable because of the extremely
large number of degrees of freedom -dof- involved. To alleviate this drawback, one possibility
lies in the use of a model reduction (based on the Karhunen-Loève decomposition —KLD—,
also known as proper orthogonal decomposition —POD—). Model reduction techniques
have been successfully applied in the finite element framework for modeling dynamic models
of distributed parameters [13]. However, in these applications several direct problems must
be solved to extract empirical functions that most efficiently represent the system. This
set of empirical eigenfunctions is used as functional basis of the Galerkin procedure to
lump the governing equation. Thus, for example, the resulting lumped parameter model
can be used to obtain the solution when the boundary conditions are changing randomly.
To avoid these preliminary costly calculations, Ryckelynck proposed in [14] to start the
resolution process from any reduced basis, using the Krylov subspaces generated by the
governing equation residual for enriching the approximation basis, at the same time that a
proper orthogonal decomposition extracts relevant information in order to maintain the low
order of the approximation basis. This procedure has been detailed in the previous section.
Moreover, there is an appropriate choice of a reduced number of weighting functions able
to solve the problem efficiently which leads to the “A priori hyper-reduction” approach
[14], but this approach, summarized previously, is not concerned in the present section.
Until now, the “a priori” model reduction techniques were proposed and applied in the
context of problems defined in a fixed mesh. In this section we analyze its application in
problems defined in domains evolving in time and two important aspects will be addressed:
(i) the effects of remeshing, with the related field projections; and (ii) the case of evolution
problems described with a set of nodes whose position evolves in time.

The plan of this part is the following: we start with a revision of the natural element discretisation technique (the finite element one is assumed to be well known by the reader) which as previously argued, allows to simulate large domains evolutions keeping the same set of nodes. In section 3 we introduce the main ideas related to the Karhunen-Loève decomposition as well as to the main ideas associated with the adaptive modeling allowing an “a priori” model reduction using a discrete formalism instead the continuous and formal description employed in the previous section. Finally, all the ideas will be illustrated considering some numerical examples involving coupled non-linear advection equations defined in domains evolving in time.

3.2 The α-Natural Element Method

In the last decade considerable research efforts have been paid to the development of a series
of novel numerical tools that have been referred as meshless or meshfree methods. These
methods do not need explicit connectivity information, as required in standard FEM. The
geometrical information is generated in a process transparent to the user, alleviating the
pre-processing stage of the method. They also present outstanding advantages in modelling complex phenomena, such as large deformation problems, forming processes, fluid flow, etc, where traditional and more experienced techniques, like the FEM, fail due to the need of excessive remeshing.

The Natural Element Method (NEM) is one of the latest meshless techniques applied in the field of linear elastostatics. It has unique features among meshless Galerkin methods, such as interpolant character of shape functions and exact application of essential boundary conditions (see the review paper [5]). These capabilities and its inherent meshless character make the NEM an appealing choice also for its application in the simulation of fluid flows. The NEM is based on the natural neighbour interpolation scheme, which in turn relies on the concepts of Voronoi diagrams and Delaunay triangulations (see Figure 2), to build Galerkin trial and test functions. These are defined as the natural neighbour coordinates (also known as Sibson’s coordinates) of the point under consideration, that is, with respect to Figure 3, the value at point $\underline{x}$ of the shape function associated with the node 1, is defined by:

$$N_1(\underline{x}) = \frac{Area(abfe)}{Area(abcd)} \qquad (40)$$

These functions are used to build the discrete system of equations arising from the application of the Galerkin method in the usual way. It has been proved, that the angles of the Delaunay triangles are not influencing the quality of the results, in opposition to the FEM. In addition, the NEM has interesting properties such as linear consistency and smoothness of the shape functions (natural neighbour coordinates are $C^1$ everywhere except at the nodes, where they are $C^0$). But perhaps the most interesting property of the Natural Element Method is the Kronecker delta property, i.e. $N_i(\underline{x}_j) = \delta_{ij}$. In opposition to the vast majority of meshless methods, the NEM shape functions are strictly interplant. This property allows an exact reproduction of linear (even bilinear in some 3D cases) displacement fields on the boundary of convex domains, since the influence of interior points vanishes along convex boundaries. This is not true in non-convex boundaries, where some specific treatment is required. The alpha-shape concept allows circumventing this difficulty when it is used in the context of a natural neighbour interpolation [4].

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Delaunay triangulation and Voronoi diagram of a cloud of points. On the right, an example of a degenerate distribution of nodes, with the two possible triangulations depicted. In this last case, four points lie in the same circumcircle and thus no single triangulation exists</figcaption></figure>

The application of the NEM for complex fluid flow simulations has been recently proposed in [12]. The main advantage of using the NEM in the framework of an updated Lagrangian formulation for simulating free or moving surface flows is the fact that the nodal position can be updated from the flow kinematics, without remeshing requirements, allowing the accurate description of large transformations and the history effects. In the same way it can be applied to solve, with great accuracy, advection-diffusion problems in
<figure><img src="image_4.png" /><figcaption>Figure 3. Definition of the Natural Neighbour coordinates of a point $\underline{x}$</figcaption></figure>

an operator splitting framework. Thus, the advection term governs the nodal position updating, whereas the diffusion equation is solved in the just updated domain, as illustrated later in this paper.

## 3.3 Reduced Order Modelling: A Discrete Approach

### 3.3.1 Introduction: the Karhunen-Loève decomposition

We assume that the evolution of a certain field $u(\underline{x}, t)$ is known (being its evolution governed by a PDE). In practical applications, this field is expressed in a discrete form, that is, it is known at the nodes of a spatial mesh and for some times $u(\underline{x}_i, t^n) \equiv u_i^n$. We can also write introducing a time discretisation $u^n(\underline{x}) \equiv u(\underline{x}, t = n\Delta t)$; $\forall n \in [1, \dots, P]$. The main idea of the Karhunen-Loève (KL) decomposition is how to obtain the most typical or characteristic structure $\phi(\underline{x})$ among these $u^n(\underline{x})$ $\forall n$. This is equivalent to obtaining a function $\phi(\underline{x})$ that maximizes $\lambda$ defined by

$$ \lambda = \frac{\sum_{n=1}^{n=P} \left[ \sum_{i=1}^{i=N} \phi(\underline{x}_i) u^n(\underline{x}_i) \right]^2}{\sum_{i=1}^{i=N} (\phi(\underline{x}_i))^2} \quad (41) $$

The maximization ($\delta\lambda = 0$) leads to:

$$ \sum_{n=1}^{n=P} \left[ \left( \sum_{i=1}^{i=N} \tilde{\phi}(\underline{x}_i) u^n(\underline{x}_i) \right) \left( \sum_{j=1}^{j=N} \phi(\underline{x}_j) u^n(\underline{x}_j) \right) \right] = \lambda \sum_{i=1}^{i=N} \tilde{\phi}(\underline{x}_i) \phi(\underline{x}_i); \quad \forall \tilde{\phi} \quad (42) $$

which can be rewritten in the form

$$ \sum_{i=1}^{i=N} \left\{ \sum_{j=1}^{j=N} \left[ \sum_{n=1}^{n=P} u^n(\underline{x}_i) u^n(\underline{x}_j) \phi(\underline{x}_j) \right] \tilde{\phi}(\underline{x}_i) \right\} = \lambda \sum_{i=1}^{i=N} \tilde{\phi}(\underline{x}_i) \phi(\underline{x}_i); \quad \forall \tilde{\phi} \quad (43) $$

Defining the vector $\underline{\phi}$ such that its $i$-component is $\phi(\underline{x}_i)$, Eq. (43) takes the following matrix form

$$ \underline{\tilde{\phi}}^T \underline{k} \underline{\phi} = \lambda \underline{\tilde{\phi}}^T \underline{\phi}; \quad \forall \underline{\tilde{\phi}} \Rightarrow \underline{k} \underline{\phi} = \lambda \underline{\phi} \quad (44) $$

where the two points correlation matrix is given by

$$ k_{ij} = \sum_{n=1}^{n=P} u^n(\underline{x}_i) u^n(\underline{x}_j) \Leftrightarrow \underline{k} = \sum_{n=1}^{n=P} \underline{u}^n(\underline{u}^n)^T \quad (45) $$
which is symmetric and positive definite. If we define the matrix $\underline{Q}$ containing the discrete field history:

$$ \underline{Q} = \begin{pmatrix} u_1^1 & u_1^2 & \cdots & u_1^n \\ u_2^1 & u_2^2 & \cdots & u_2^n \\ \vdots & \vdots & \ddots & \vdots \\ u_N^1 & u_N^2 & \cdots & u_N^n \end{pmatrix} \qquad (46) $$

is direct to verify that the matrix $\underline{k}$ in Eq. (44) results

$$ \underline{k} = \underline{Q} \underline{Q}^T \qquad (47) $$

where the diagonal components are given by

$$ k_{ii} = (\underline{Q} \underline{Q}^T)_{ii} = \sum_{j=1}^{j=P} (u_i^j)^2 \qquad (48) $$

Thus, the functions defining the most characteristic structure of $u^n(\underline{x})$ are the eigenfunctions $\phi_k(\underline{x}) \equiv \underline{\phi}_k$ associated with the highest eigenvalues.

### 3.3.2 A posteriori reduced modelling

If some direct simulations have been carried out, we can determine $u(\underline{x}_i, t^n) \equiv u_i^n$, $\forall i \in [1, \dots, N] \ \forall n \in [1, \dots, P]$, and from these the $r$ eigenvectors related to the $r$-highest eigenvalues $\underline{\phi}_k = \phi_k(\underline{x}_i)$, $\forall i \in [1, \dots, N]$, $\forall k \in [1, \dots, r]$ (with $r \ll N$). Now, we can try to use these $r$ eigenfunctions for approximating the solution of a problem slightly different to the one that has served to define $u(\underline{x}_i, t^n)$. For this purpose we need to define the matrix $\underline{A}$

$$ \underline{A} = \begin{pmatrix} \phi_1(\underline{x}_1) & \phi_2(\underline{x}_1) & \cdots & \phi_r(\underline{x}_1) \\ \phi_1(\underline{x}_2) & \phi_2(\underline{x}_2) & \cdots & \phi_r(\underline{x}_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_1(\underline{x}_N) & \phi_2(\underline{x}_N) & \cdots & \phi_r(\underline{x}_N) \end{pmatrix} \qquad (49) $$

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

$$ \underline{H} \underline{U} = \underline{F} \qquad (50) $$

Obviously, in the case of evolution problems $\underline{F}$ contains the contribution of the solution at the previous time step.

Then, assuming that the unknown vector contains the nodal degrees of freedom, it can be expressed as

$$ \underline{U} = \sum_{i=1}^{i=r} a_i \underline{\phi}_i = \underline{A} \underline{a} \qquad (51) $$

it results

$$ \underline{H} \underline{U} = \underline{F} \Rightarrow \underline{H} \underline{A} \underline{a} = \underline{F} \qquad (52) $$

and multiplying both terms by $\underline{A}^T$ it results

$$ \underline{A}^T \underline{H} \underline{A} \underline{a} = \underline{A}^T \underline{F} \qquad (53) $$

which proves that the final system of equations is of low order, i.e. the dimensions of $\underline{A}^T \underline{H} \underline{A}$ are $r \times r$, with $r \ll N$, and the dimensions of both $\underline{a}$ and $\underline{A}^T \underline{F}$ are $r \times 1$.
*Remark 3.1.* Eq.(53) can be also derived introducing the approximation (51) into the PDE Galerkin form.

*Remark 3.2.* Ryckelynck proposed in [14] another reduction in the number of weighting functions required to evaluate accurately equation (53). He called this new technique an “a priori” hyper-reduction strategy. Even if this strategy could be easily applied in the examples that we present later, we prefer, for the sake of simplicity, avoid its introduction in those examples because it has been discussed deeply in the first part of this work.

### 3.3.3 Adaptivity via an a priori model reduction: discrete approach

In order to compute reduced model solutions without an a priori knowledge, Ryckelynck proposed in [14] 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. Now, we compute S iterations of the evolution problem using the reduced model (53) without changing the approximation basis $A^{(0)}$ (the superscript indicates that this is the first approximation basis used). For a more general description we consider that the approximation basis has been updated m times until now, that is, the approximation basis at the present time is $A^{(m)}$. Now, we compute S steps of the reduced model without changing the approximation basis. After each S iterations, the complete discrete system (52) is constructed, and the residual $\underline{R}$ evaluated:

$$ \underline{R} = \underline{H} \underline{U} - \underline{F} = \underline{H} \underline{A}^{(m)} \underline{a}^{(m)} - \underline{F} \quad (54) $$

If the norm of the residual verifies $||\underline{R}|| < \varepsilon$, with $\varepsilon$ a threshold value small enough, we can continue for other S iterations using the same approximation basis $A^{(m)}$. On the contrary, if the residual norm is too large, $||\underline{R}|| > \varepsilon$, we put $S = S/2$ and recompute. When the residual becomes small enough, i.e. $||\underline{R}|| < \varepsilon$ the reduced approximation basis is enriched using some Krylov's subspaces $\{\underline{R}, \underline{H}, \underline{R}, \underline{H}^2 \underline{R}, \dots\}$. At this time the reduced solution $\underline{a}_n^{(m)}$ is stocked ($n$ represents the corresponding time step).

Now, we compute the most representative information extracted from the reduced solutions $\underline{a}_n^{(m)}$ previously stocked. The superscript (m) indicates that these reduced order solutions are expressed in the basis $A^{(m)}$. Now, applying the Karhunen-Loève decomposition to the solution evolution represented for these vectors $\underline{a}_n^{(m)}$ ($\forall n$) we obtain the most representative eigenvectors defining the matrix $\underline{V}$.

Then the evolution process can continue for other S time steps, using the enriched basis defined by: $A^{(m+1)} = \{A^{(m)}\underline{V}, \underline{R}, \underline{H}, \underline{R}, \underline{H}^2\underline{R}\}$ (in our simulation we consider only the first three Krylov's subspaces).

After each reduced basis modification, the previous reduced solutions that have been stocked $\underline{a}_n^{(m)}$ must be projected into the new basis. Thus, we can write:

$$ \underline{a}_n^{(m+1)} = [(\underline{A}^{(m+1)})^T \underline{A}^{(m+1)}]^{-1} (\underline{A}^{(m+1)})^T \underline{A}^{(m)} \underline{a}_n^{(m)}, \quad \forall n \quad (55) $$

*Remark 3.3.* Eq.(55) is formally equivalent to Eq.(20).

*Remark 3.4.* The application of the Karhunen-Loève decomposition to the reduced order solutions instead to the nodal description has two advantages: (i) the eigenvalue problem has a lower dimension, and (ii) as the functions in $A$ verify the problem boundary conditions, then any low order solution $\underline{a}$ determines a nodal solution $\underline{A} \underline{a}$ verifying these boundary conditions.
Remark 3.5. In the numerical examples that follow, we will consider only one time interval divided in S steps. Thus, with respect to the procedure just described, we will consider only the approximation basis enrichment, to conclude about the convergence of this technique in the framework of non-linear coupled advection-diffusion equations.

## 3.4 Numerical Examples

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

$$ \underline{v} = \begin{pmatrix} u(x, y) \\ v(x, y) \end{pmatrix} = \begin{pmatrix} \dot{\gamma}x \\ -\dot{\gamma}y \end{pmatrix} \qquad (56) $$

being its length and width at time t, $2L_t = 2(L_0e^{\gamma t})$ and $2H_t = 2(H_0e^{-\gamma t})$ respectively, being $\Omega_t = [-L_t, L_t] \times [-H_t, H_t]$. Thus it results $|\Omega_t| = |\Omega_0|$.

We consider the following non-linear and coupled advection-diffusion problem defined in this domain:

$$ \left\{ \begin{aligned} \frac{dT}{dt} &= \alpha \Delta T - \beta \gamma(T) C \\ \frac{dC}{dt} &= D \Delta C - \gamma(T) C \end{aligned} \right. \qquad (57) $$

Remark 3.6. This problem can be viewed as modeling the heat transfer in presence of an exothermic reaction whose kinetics depends itself on the concentration and on the temperature in a non-linear way.

Now, we can proceed to discretize the problem defined by Eq. (57) using an splitting operator technique as well as an implicit Galerkin finite or natural element technique (the problematic related to the mesh distortion will be discussed later) for solving the resulting diffusion problems in the updated geometry. Firstly, we will consider a fixed point linearization. Thus, Eq.(57) can be solved in two steps:

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

$$ \left\{ \begin{aligned} \underline{x}_i^{n+1} &= \underline{x}_i^n + \underline{v}_i^n \Delta t \\ T_i^{n+\frac{1}{2}}(\underline{x}_i^{n+1}) &= T_i^n(\underline{x}_i^n) \\ C_i^{n+\frac{1}{2}}(\underline{x}_i^{n+1}) &= C_i^n(\underline{x}_i^n) \end{aligned} \right. \qquad (58) $$

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

$$ \begin{aligned} & (\underline{T}^{*T} \underline{C}^{*T}) \left(\frac{\underline{M}}{\underline{0}}, \frac{\underline{0}}{\underline{M}}\right) \left(\frac{\underline{T}^{n+1}}{\underline{C}^{n+1}}\right) = (\underline{T}^{*T} \underline{C}^{*T}) \left(\frac{\underline{M}}{\underline{0}}, \frac{\underline{0}}{\underline{M}}\right) \left(\frac{\underline{T}^{n+\frac{1}{2}}}{\underline{C}^{n+\frac{1}{2}}}\right) - \\ & - \Delta t (\underline{T}^{*T} \underline{C}^{*T}) \left(\frac{\underline{\alpha K}}{\underline{0}}, \frac{\underline{0}}{\underline{D K}}\right) \left(\frac{\underline{T}^{n+1}}{\underline{C}^{n+1}}\right) - \Delta t (\underline{T}^{*T} \underline{C}^{*T}) \left(\frac{\underline{0}}{\underline{0}}, \frac{\underline{\beta G}}{\underline{G}}\right) \left(\frac{\underline{T}^{n+1}}{\underline{C}^{n+1}}\right) \end{aligned} \qquad (59) $$

where the particular form of matrix $\underline{G}$ depends on the linearisation of the non linear term $\gamma(T)$.
### 3.4.1 Uncoupled model analysis

We start our analysis considering the uncoupled model resulting when in Eq.(59) one assumes $\gamma(T) = 0$. Now, 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 \underline{n}|_{\partial\Omega_t} = 0$ (being $\underline{n}$ the unit outwards vector defined on the domain boundary), and the following initial temperature is considered in $\Omega_0$:

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

The numerical simulation is carried out with the choice of parameters grouped in Table 1.

<table><caption>**Table 1.** Simulation parameters considered in the numerical examples</caption><thead><tr><th>N</th><th>t<sub>max</sub></th><th>L = 0</th><th>T<sub>0</sub></th><th>Δt</th><th>˙γ</th><th>α</th><th>β</th><th>D</th><th>γ</th><th>C<sub>0</sub></th></tr></thead><tbody><tr><td>101 × 10</td><td>5</td><td>π</td><td>2</td><td>0.1</td><td>0.1</td><td>0.5</td><td>1</td><td>1</td><td>0</td><td>2</td></tr></tbody></table>

*Remark 3.7.* The number of nodes assumed in the y-direction is irrelevant due to the unidirectional heat transfer resulting from the initial temperature field (60). This fact facilitates the model validation as well as the solution representation.

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

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

$$\underline{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} \qquad (61)$$

It is very important to notice that these vectors are associated with different nodal distributions. Now, we solve the eigenvalue problem (44)

$$(\underline{Q} \underline{Q}^T) \underline{\phi} = \lambda \underline{\phi} \qquad (62)$$

from which $N$ eigenvalues result, each one related to one eigenvector. As previously described, we select only the eigenvectors related to eigenvalues in the interval defined by the highest eigenvalue and $10^{-8}$ times this value, i.e. $\{\underline{\phi}_1, \dots, \underline{\phi}_r\}$ such that $\lambda_k \in [10^{-8}\lambda_1, \lambda_1]$, $\forall k \in [1, \dots, r]$, being $\lambda_1$ the highest eigenvalue. In the problem just simulated there are only two eigenvalues in this interval, and in consequence the whole time evolution could be represented, very accurately, as a linear combination of the two associated eigenvectors that we denote by $\underline{\phi}_1$ and $\underline{\phi}_2$. We emphasize the fact that both vectors are associated to the nodes but not to a particular nodal distribution. Moreover, in order to represent the initial condition we propose to add to both eigenvectors, the vector form of the initial condition that we denote by $\underline{\phi}_0 = \underline{T}^0$. Thus the low order approximation basis
<figure><img src="image_5.png" /><figcaption>Figure 4. Initial temperature distribution (stars) and temperature profiles (along the x-axis) at times 5Δt, 10Δt, ..., 50Δt</figcaption></figure>

which consist of the three functions {$\phi_0, \phi_1, \phi_2$} are no more orthogonal. Figure 5 depicts these functions after normalization that we represent on the final geometry. Now, we are going to prove that these three approximation functions are enough to represent accurately the whole temperature evolution. For this purpose we write:

$$ \underline{T} = \sum_{i=0}^{i=2} a_i(t) \underline{\phi}_i = \underline{A} \underline{a} = \begin{pmatrix} \phi_0(\underline{x}_1) & \phi_1(\underline{x}_1) & \phi_2(\underline{x}_1) \\ \vdots & \vdots & \vdots \\ \phi_0(\underline{x}_N) & \phi_1(\underline{x}_N) & \phi_2(\underline{x}_N) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \\ a_2 \end{pmatrix} \quad (63) $$

that we introduce in the first relation of Eq. (59), which involves the temperature field, to perform the following reduced order linear system:

$$ \underline{T}^{*T}[\underline{M} + \alpha \Delta t \underline{K}] \underline{T}^{n+1} = \underline{T}^{*T} \underline{M} \underline{T}^n \Rightarrow \underline{a}^{*T} \underline{A}^T [\underline{M} + \alpha \Delta t \underline{K}] \underline{A} \underline{a}^{n+1} = \underline{a}^{*T} \underline{A}^T \underline{M} \underline{A} \underline{a}^n \quad (64) $$

after which the reduced order solution $\underline{a}^{n+1}$ can be computed from the previous one according to

$$ \underline{a}^{n+1} = [\underline{A}^T[\underline{M} + \alpha \Delta t \underline{K}]\underline{A}]^{-1} \underline{A}^T \underline{M} \underline{A} \underline{a}^n \quad (65) $$

where the initial condition is expressed in the reduced basis approximation by

$$ (\underline{a}^0)^T = (1 \ 0 \ 0) \quad (66) $$

In Figure 6 we depict the computed solution using the reduced approximation basis versus the one computed using the natural element method. As we can notice the evolution of the solution is perfectly described despite of the few number of functions involved in the reduced order approximation.
<figure><img src="image_6.png" /><figcaption>Figure 5. Normalized approximation functions: the continuous line represents the initial condition, and the dashed and dotted curves the two problem eigenvectors</figcaption></figure>

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, i.e.:

$$ \underline{T} = \underline{\phi}_0(\underline{x})a(t) = \begin{pmatrix} \phi_0(\underline{x}_1) \\ \vdots \\ \phi_0(\underline{x}_N) \end{pmatrix} a(t) = \underline{A}^{(0)}\underline{a} \qquad (67) $$

Introducing the expression of the approximation basis $\underline{A}^{(0)}$ given by Eq.(67) in the equation governing the evolution of the reduced order solution (Eq.(65)), and taking into account that the verification of the initial condition implies that $a(t = 0) = 1$, we can compute its evolution, i.e. $a(t)$, and from this the evolution of $T(\underline{x}, t)$ ($T(\underline{x}, t) = \underline{A}^{(0)}\underline{a}$). Figure 7 compares the reduced order solutions computed at different times (stars curves) with the reference ones depicted in Figure 4 (continuous curves).

Now, from the final reduced order solution at time $50\Delta t$ (Figure 7(d), stars curve), which we denote by $\underline{T}^{(0)}(t = 50\Delta t)$, we can compute the residual that allows us to enrich the approximation basis using the Krylov's subspaces. Thus, the residual derives automatically from Eq.(59)

$$ \underline{R} = [\underline{M} + \alpha\Delta t\underline{K}]\underline{T}^{(0)}(t = 50\Delta t) - \underline{M}\underline{T}^{(0)}(t = 49\Delta t) \qquad (68) $$

which allows to compute the successive Krylov's subspaces. Thus, if we define the matrix $\underline{H}$ by:

$$ \underline{H} = \underline{M}^{-1}[\underline{M} + \alpha\Delta t\underline{K}] \qquad (69) $$
<figure><img src="image_7.png" /><figcaption>Figure 6. Reduced order solutions at different times —5Δt, 10Δt, ..., 50Δt— computed using the reduced order approximation basis (cross symbol) versus the reference ones computed using an approximation based on the whole set of nodes</figcaption></figure>

then, the different Krylov's subspaces $\underline{KR}_k$ can be written as

$$ \underline{KR}_k = \underline{H}^k \underline{R}, \forall k \ge 0 \qquad (70) $$

The residual in our problem (which defines the first Krylov's subspace) is depicted in Figure 8 (in fact the successive Krylov's subspaces are very close to the first one). 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 at each time (that is, in each domain).

As the residual norm is higher than a small enough value $\varepsilon$ ($\varepsilon = 10^{-6}$) we must re-computed the whole evolution using the enriched basis. Thus, we consider the reduced approximation enriched basis $\underline{A}^{(1)}$, which is obtained, adding the residual just computed to $\underline{A}^{(0)}$. Thus the new approximation results:

$$ \underline{T} = \begin{pmatrix} \phi_0(\underline{x}_1) & R(\underline{x}_1) \\ \vdots & \vdots \\ \phi_0(\underline{x}_N) & R(\underline{x}_N) \end{pmatrix} \begin{pmatrix} a_0 \\ a_1 \end{pmatrix} = \underline{A}^1 \underline{a} \qquad (71) $$

that introduced as previously in Eq.(65) allows to recompute the evolution of the solution in $[0, t_{max} = 50\Delta t]$ (with $\underline{a}^T(t = 0) = (1 \ 0)$), and the new residual at the last time step according to Eq.(68) substituting $\underline{T}^{(0)}(t = 50\Delta t)$ by the current $\underline{T}^{(1)}(t = 50\Delta t)$. In our case the evolution obtained is very close to the one depicted in Figure 6, and the norm of the new residual is $\|R\| = 10^{-10}$ which proves that with only two approximation functions the evolution can be represented with very high accuracy. Moreover we have
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Reduced order solutions computed using the approximation basis consisting on the initial condition (stars) versus de reference ones (continuous curve), 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>

proved that the enrichment technique based on the Krylov's subspaces, originally proposed in the Ryckelynck works, is able to adapt also the approximation basis in few iterations in problems involving large domain transformations.

### 3.4.2 Projection induced by remeshing

When one proceeds in the framework of the natural element discretisation, large domain transformations can be accounted without any remeshing requirement, instead the frequent remeshing required when one proceeds using a finite element strategy. However, even in the natural element case, some times we need to add new nodes during the problem resolution, according to a suitable error estimator or indicator, to improve the solution representation, and consequently the computed solution accuracy. In the same way, some times, one could wish reduce de number of nodes, because even if there is not a direct incidence in the size of the linear problem when a model reduction is applied, they imply an integration necessity in the discretisation of the variational formulation. One possibility for alleviating this arduous task lies in the use of a selective integration as proposed in Ryckelynck [14] and summarized in the first part of this paper. In any case we will analyze in this section the incidence of a remeshing in the low order approximation basis, and we prove that if a
<figure><img src="image_9.png" /><figcaption>Figure 8. Computed residual at the last time step related to the low order basis $\underline{A}^{(0)}$</figcaption></figure>

basis projection follows the remeshing operation, the results quality remains unchanged.

To illustrate this procedure, we consider that the evolution is started at time $t = 0$ using the previous reduced approximation basis $\underline{A}^{(1)}$, from the initial condition $\underline{a}^T(t = 0) = (1 \ 0)$. At time $t = 25\Delta t$ the nodes are repositioned according to the expression (72) (other possibilities have been also tested):

$$x_i^* = x_i + (L_t - x_i) \frac{|x_i|}{2L_t}, \quad \forall i \in [1, \dots, N] \qquad (72)$$

without changes in the y-coordinate.

In fact, if despite of the remeshing operation the approximation basis is not updated, the approximation functions expressed from the new mesh result significantly different to the initial ones (Figure 9(a)), and consequently, large deviations appear in the final solution with respect to the reference one, as noticed in Figure 9(b).

To avoid the deviations noticed in Figure 9(b), we propose to apply a projection of the approximation functions between the old and new meshes. Several possibilities exist to perform this projection, but here we consider the simplest one, that is defined by the following relation:

$$\begin{pmatrix} \phi_j^*(\underline{x}_1^*) \\ \phi_j^*(\underline{x}_2^*) \\ \vdots \\ \phi_j^*(\underline{x}_N^*) \end{pmatrix} = \begin{pmatrix} N_1(\underline{x}_1^*) & N_2(\underline{x}_1^*) & \cdots & N_N(\underline{x}_1^*) \\ N_1(\underline{x}_2^*) & N_2(\underline{x}_2^*) & \cdots & N_N(\underline{x}_2^*) \\ \vdots & \vdots & \ddots & \vdots \\ N_1(\underline{x}_N^*) & N_2(\underline{x}_N^*) & \cdots & N_N(\underline{x}_N^*) \end{pmatrix} \begin{pmatrix} \phi_j(\underline{x}_1) \\ \phi_j(\underline{x}_2) \\ \vdots \\ \phi_j(\underline{x}_N) \end{pmatrix}; \quad j = 1, 2 \qquad (73)$$

where $N_k$ are the shape functions related to the old nodal distribution. The previous system
<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) Representation of the normalized approximation functions before (continuous) and after the nodes repositioning (dashed); (b) final computed field at $t = 50\Delta t$ when a remeshing is applied at time $t = 25\Delta t$ without changes in the reduced approximation basis: computed solution (dashed) compared to the reference one (continuous line)</figcaption></figure>

can be expressed in the compact form:

$$ \underline{J} \underline{\phi}_j = \underline{\phi}_j^*; \quad j = 1, 2. \qquad (74) $$

Thus, after each remeshing the change of basis matrix $\underline{J}$ must be computed, and the approximation basis projected in the new one according to:

$$ \underline{A}^* = \underline{J} \underline{A} \qquad (75) $$

that replaces the expression of $\underline{A}$ in Eq. (65) up to the remeshing step, until the new remeshing which involves a new projection.

If we proceed in this form, then we obtain the results shown in Figure 10 instead the previously obtained (without the projection) depicted in Figure 9.

### 3.5 An Example of Non-Linear Coupled Models

In this section we consider the coupled formulation given by Eq. (59) when the parameter $\delta(T)$ is no more null, and a fixed-point linearization is assumed for $\gamma(T) = \gamma T$. Thus, for the simulation we consider the parameters grouped in Table 1 with $\gamma = 1$. Moreover, the initial concentration field is assumed in the form

$$ C(x, y, t = 0) = C_0(1 + \cos(x)) \qquad (76) $$

We assume a reduced coupled approximation basis such that

$$ \left[ \frac{T^{n+1}}{C^{n+1}} \right] = \hat{A} \hat{a} \qquad (77) $$

which introduced in Eq. (59) allows the coupled solution updating:

$$ \hat{a}^{n+1} = \left[ \hat{A}^T [\hat{M} + \Delta t \hat{K} + \Delta t \hat{G}^n] \hat{A} \right]^{-1} \hat{A}^T \hat{M} \hat{A} \hat{a}^n \qquad (78) $$
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) Representation of the non-normalized approximation functions before (stars) and after the nodes repositioning (continuous); (b) final computed field at $t = 50\Delta t$ when a remeshing is applied at time $t = 25\Delta t$ and then the associated projection of the reduced approximation basis: computed solution (stars) compared to the reference one (continuous)</figcaption></figure>

where

$$ \hat{M} = \begin{pmatrix} \underline{M} & \underline{0} \\ \underline{0} & \underline{M} \end{pmatrix} \qquad (79) $$

$$ \hat{K} = \begin{pmatrix} \alpha \underline{K} & \underline{0} \\ \underline{0} & \underline{D K} \end{pmatrix} \qquad (80) $$

$$ \hat{G}^n = \begin{pmatrix} \underline{0} & \beta \underline{G}^n \\ \underline{0} & \underline{G}^n \end{pmatrix} \qquad (81) $$

*Remark 3.8.* In the case of $\gamma(T) = \gamma = 1$ the problem becomes linear and the matrix $\underline{G}^n$ is reduced to the mass matrix, i.e. $\underline{G}^n = \underline{M}$.

Now, we start the evolution process with the following reduced coupled approximation basis $\hat{A}^{(0)}$ and reduced order coupled initial condition $\hat{a}(t = 0)$:

$$ \hat{A}^{(0)} = \begin{pmatrix} \underline{T}^0 & \underline{0} \\ \underline{0} & \underline{C}^0 \end{pmatrix} \qquad (82) $$

where $\underline{T}^0$ and $\underline{C}^0$ are the vector representation of both nodal temperatures and concentrations. Thus, the natural election of $\hat{a}(t = 0)$ results

$$ \hat{a}(t = 0) = \begin{pmatrix} 1 \\ 1 \end{pmatrix} \qquad (83) $$

The temperature and concentration reduced order solution are depicted in Figure 11 at times $5\Delta t$, $20\Delta t$, $35\Delta t$ and $50\Delta t$.

In Figure 12 we depict the first Krylov's subspace (residuals) for both the temperature (dashed curve) and the concentration (continuous curve).
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Reduced order temperatures (continuous) and concentrations (dashed) computed using the approximation basis consisting only the initial conditions, at times $5\Delta t$, $20\Delta t$, $35\Delta t$ and $50\Delta t$ (up-left, up-right, down-left and down-right respectively)</figcaption></figure>

The basis enrichment is performed by adding to $\hat{\mathcal{A}}^{(0)}$ the residual just computed to generate the next approximation basis $\hat{\mathcal{A}}^{(1)}$, which will be used to recompute the whole evolution process. The other Krylov's subspaces are not introduced because they result very close to the first one. This adaptive procedure, in which at the end of the evolution process the approximation basis is updated using the Krylov's subspaces continues until the convergence, that is, until obtaining a small enough residual norm. Figure 13 illustrates the evolution solution after reaching convergence in the adaptive procedure (this convergence is achieved in some - two or three - iterations despite the use of a fixed-point linearization of the non-linear term). This solution is very accurate with respect to the reference one computed by solving the problem using the finite or the natural element representation. We can notice in Figure 13 the non-linear coupling effect which implies higher reaction rates as the temperature increases. Thus, in the zones with higher temperature, we observe a faster evolution of the concentration, which reaches its lowest values.

Other possibility for treating the non-linearity lies in the use of a Newton technique. In this case at each time step the solution is searched from the one computed at the previous iteration using the iteration algorithm proposed in the next paragraphs.
<figure><img src="image_13.png" /><figcaption>Figure 12. Computed residuals at the last time step related to the low order basis $\hat{\underline{A}}^{(0)}$ for both temperatures (dashed) and concentrations (continuous)</figcaption></figure>

We start from Eq. (59) that we write in a compact form

$$[\underline{\hat{M}} + \Delta t \underline{\hat{K}} + \Delta t \underline{\hat{G}}(\underline{TC}^{n+1})] \underline{TC}^{n+1} - \underline{\hat{M}} \underline{TC}^{n} = 0 \quad (84)$$

where the $TC$ is the vector containing the nodal temperatures and concentrations.

Now assuming that $TC^n$ and $TC^{n+1}$ will be expressed using the reduced approximation basis, we can also consider that if we are looking for the value of $TC^{n+1}$ starting from an estimation (that could consist in the solution at the previous time step) by adding the pertinent corrections, then that correction can be expressed using the same reduced order approximation basis. Thus, finally, we can write:

$$[\underline{\hat{A}}^T [\underline{\hat{M}} + \Delta t \underline{\hat{K}} + \Delta t \underline{\hat{G}}^{n+1(s)}] \underline{\hat{A}}] \underline{\Delta \hat{a}}^{n+1(s)} = -[\underline{\hat{A}}^T [\underline{\hat{M}} + \Delta t \underline{\hat{K}} + \Delta t \underline{\hat{G}}^{n+1(s)}] \underline{\hat{A}}] \underline{\hat{a}}^{n+1(s)} + \underline{\hat{A}}^T \underline{\hat{M}} \underline{\hat{A}} \underline{\hat{a}}^n \quad (85)$$

with

$$\underline{\hat{a}}^{n+1(s+1)} = \underline{\hat{a}}^{n+1(s)} + \underline{\Delta \hat{a}}^{n+1(s)} \quad (86)$$

and being $\underline{\hat{G}}^{n+1(s)}$ the tangent matrix related to the linearization of the non-linear term (see Eq. (57)), whose variational expression results:

$$\gamma \int C^*TC \Rightarrow \gamma \int C^*\Delta TC + \gamma \int C^*T\Delta C \quad (87)$$
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Reduced order temperatures (dashed) and concentrations (continuous) computed using the approximation basis obtained after the convergence of the adaptive enrichment procedure, at times $5\Delta t$, $20\Delta t$, $35\Delta t$ and $50\Delta t$ (up-left, up-right, down-left and down-right respectively)</figcaption></figure>

We have noticed that the non-linear algorithm, that must be solved at each time step due to the full-implicit character of the scheme, converges in two or three iterations.

# 4 APPLICATION OF MODEL REDUCTION IN THE SIMULATION OF KINETIC THEORY MODELS

## 4.1 Mechanical Modelling of Short Fiber Suspensions

Numerical modelling of non-Newtonian flows usually involves the coupling between equations of motion, which define an elliptic problem, and the fluid constitutive equation, which introduces an advection problem related to the fluid history. In short fiber suspensions (SFS) models, the extra-stress tensor depends on the fiber orientation whose evolution can be modelled from a transport problem. In all cases the flow kinematics and the fiber orientation are coupled: the kinematics of the flow governs the fiber orientation, and the presence and orientation of the fibers modify the flow kinematics.

If one uses SFS flows in material forming processes, the final fiber orientation state depends on the process and exhibits flow-induced anisotropy. Thus, we need to compute
the fiber orientation in order to predict the final mechanical properties of the composite parts, which depends strongly on the fiber orientation. Moreover, the numerical simulation of such flows becomes interesting if one want to identify their rheological parameters using some rheometric devices and an appropriate inverse technique.

This section focuses in the determination of fiber orientation assuming a known kinematics. For this purpose an accurate approach deal with a microscopic simulation using an orientation distribution function related to each material point. That formalism, called Fokker-Planck formalism, is a commonly used description of kinetic theory problems, for describing the evolution of the configuration distribution function. This function represents the probability of finding the microstructural element in a particular configuration.

In the case of a short fiber suspension, the configuration distribution function (also known as orientation distribution function) gives the probability of finding the fiber in a given direction. Obviously, this function depends on the physical coordinates (space and time) as well as on the configuration coordinates, that taking into account the rigid character of the fibers, are defined on the surface of the unit sphere. Thus, we can write $\Psi(\underline{x}, t, \underline{p})$, where $\underline{x}$ defines the position of the fiber center of mass, $t$ the time and $\underline{p}$ the unit vector defining the fiber orientation. The evolution of the distribution function is given by the Fokker-Planck equation

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial \underline{p}}(\Psi \dot{\underline{p}}) + \frac{\partial}{\partial \underline{p}}\left(D_r \frac{\partial \Psi}{\partial \underline{p}}\right) \quad (88) $$

where $d/dt$ represents the material derivative, $D_r$ is a diffusion coefficient and $\dot{\underline{p}}$ is the fiber rotation velocity. The orientation distribution function must verify the normality condition:

$$ \oint \Psi(\underline{p}) d\underline{p} = 1 \quad (89) $$

When the fibers are assumed ellipsoidal and when the suspension is dilute enough, the rotation velocity can be obtained from the Jeffery's equation

$$ \dot{\underline{p}} = \underline{\Omega} \underline{p} + k \underline{D} \underline{p} - k(\underline{p}^T \underline{D} \underline{p}) \underline{p} \quad (90) $$

where $\underline{\Omega}$ and $\underline{D}$ are the vorticity and the strain rate tensors respectively, associated with the fluid flow undisturbed by the presence of the fiber, and $k$ is a scalar which depends on the fiber aspect ratio $\lambda$ (ratio between the fiber length and the fiber diameter)

$$ k = \frac{\lambda^2 - 1}{\lambda^2 + 1} \quad (91) $$

In a former paper [3] the discretisation of the advection dominated Fokker-Planck equation, governing the fiber orientation in short fiber suspension flows, was carried out using a particle technique, where the diffusion term was modelled from random motions. It was pointed out that the number of fibers required in this stochastic simulation to describe the fiber distribution increases significantly with the diffusion coefficient $D_r$. Thus, it was argued that for practical applications the use of the particle method in the framework of a stochastic simulation, is restricted to very slight diffusion effects. When the diffusion becomes dominant, continuous approximations using a fixed mesh seem to be suitable, but in this case accurate stabilizations are required for dealing with small diffusion effects, and a lack of accuracy is noticed in the treatment of the advection dominated case.

In order to reduce degrees-of-freedom (DOF) needed for computing the orientation distribution, another particle approach [1] was used in order to solve the Fokker-Planck equation.
equation using smooth particles. The main idea of this deterministic approach lies in the introduction of the Fokker-Planck diffusion term into the advection one, which allows to proceed in a Lagrangian deterministic manner without a mesh support requirement. However, the Fokker-Planck equation is defined in a multidimensional space involving the physical and the conformation coordinates. This fact induced the necessity of using a extremely large number of particles, whit an unfavorable incidence in the method efficiency.

In this section an “a priori” model reduction technique will be introduced for reducing the number of degrees of freedom involved in the simulation.

## 4.2 Spatial Discretisation

An illustration of the technique will be presented in the 2D case where $p^T = (\cos \varphi, \sin \varphi)$. Thus the Fokker-Planck equation writes:

$$ \frac{d\Psi}{dt} = -\frac{\partial}{\partial\varphi}(\Psi\dot{\varphi}) + \frac{\partial}{\partial\varphi}\left(D_r\frac{\partial\Psi}{\partial\varphi}\right) \quad (92) $$

where $\varphi$ is the orientation angle, $\Psi$ and $\dot{\varphi}$ are functions of $(x, t, \varphi)$, and the diffusion factor $D_r$ is assumed constant. For the seek of simplicity, time will be considered as an implicit variable and material position $x$ will be omitted knowing that a distribution function will be associated to each material point. Obviously, in homogeneous field, only one distribution function will be computed. Equation (92) can be written as

$$ \frac{d\Psi}{dt} + E_0(\varphi)\Psi + E_1(\varphi)\frac{\partial\Psi}{\partial\varphi} + E_2(\varphi)\frac{\partial^2\Psi}{\partial\varphi^2} = 0 \quad (93) $$

The velocity vector has two components $v^T = (u, v)$ and $(x, y)$ are the space coordinates. Functions $E_0$, $E_1$ and $E_2$ can be deduced from the velocity gradient

$$ \mathrm{Grad}v = \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{\partial u}{\partial y} \\ \frac{\partial v}{\partial x} & \frac{\partial v}{\partial y} \end{pmatrix} \quad (94) $$

$$ \left\{ \begin{aligned} E_0(\varphi) &= -k \left( \left(\frac{\partial u}{\partial x} - \frac{\partial v}{\partial y}\right) \cos(2\varphi) + \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \sin(2\varphi) \right) \\ E_1(\varphi) &= \frac{1}{2} \left( \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} + k \left[ \left(\frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}\right) \cos(2\varphi) + \left(\frac{\partial v}{\partial y} - \frac{\partial u}{\partial x}\right) \sin(2\varphi) \right] \right) \\ E_2(\varphi) &= E_2 = -D_r \end{aligned} \right. \quad (95) $$

Firstly, the problem is formulated in the Finite Element framework using a weighting function $\Psi^*$.

$$ \int_{-\pi}^{\pi} \frac{d\Psi}{dt} \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_0(\varphi) \cdot \Psi \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_1(\varphi) \cdot \frac{\partial\Psi}{\partial\varphi} \cdot \Psi^* d\varphi \\ + \int_{-\pi}^{\pi} E_2(\varphi) \cdot \frac{\partial^2\Psi}{\partial\varphi^2} \cdot \Psi^* d\varphi = 0 \quad (96) $$

Using a linear interpolation in each 2-nodes 1D-element, we can write:

$$ \Psi^{e*}(\varphi) = \sum_{i=1}^{2} \Psi_{i}^{*}(\varphi) = \sum_{i=1}^{2} N_{i}(\varphi) \Psi_{i}^{e*} \quad (97) $$
Where $\Psi_i^e$ and $\Psi_i^{e*}$ are the values at node $i$ of $\Psi$ and $\Psi^*$ respectively, and $N_i(\varphi)$ is the associated linear shape function which takes a unit value at the node $i$ vanishing at the other nodal positions.

Assuming the solution periodicity Eq.(96) can be rewritten after integration by parts, as:

$$ \int_{-\pi}^{\pi} \frac{d\Psi}{dt} \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_0(\varphi) \cdot \Psi \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_1(\varphi) \cdot \frac{\partial\Psi}{\partial\varphi} \cdot \Psi^* d\varphi - \int_{-\pi}^{\pi} E_2(\varphi) \cdot \frac{\partial\Psi}{\partial\varphi} \cdot \frac{\partial\Psi^*}{\partial\varphi} d\varphi = 0 \quad (98) $$

Eq.(98) exhibits an advection-diffusion character. The Finite Element Method needs a special treatment to avoid numerical instabilities induced by the convection term. A non-consistent upwinding formulation is considered here, which modifies the weighting function related to the advection term $\bar{\Psi}^*$ as described later. Thus, the stabilized variational formulation results:

$$ \int_{-\pi}^{\pi} \frac{d\Psi}{dt} \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_0(\varphi) \cdot \Psi \cdot \Psi^* d\varphi + \int_{-\pi}^{\pi} E_1(\varphi) \cdot \frac{\partial\Psi}{\partial\varphi} \cdot \bar{\Psi}^* d\varphi - \int_{-\pi}^{\pi} E_2(\varphi) \cdot \frac{\partial\Psi}{\partial\varphi} \cdot \frac{\partial\Psi^*}{\partial\varphi} d\varphi = 0 \quad (99) $$

where in each element

$$ \bar{\Psi}^{e*}(\varphi) = \sum_{i=1}^{2} \bar{\Psi}_{i}^{*}(\varphi) \quad (100) $$

with

$$ \bar{\Psi}_i^*(\varphi) = N_i(\varphi) + \frac{\beta h}{2} \frac{\partial N_i(\varphi)}{\partial \varphi} \quad (101) $$

where $h$ is the element size and $\beta$ is the upwinding parameter whose optimal value results:

$$ \beta = \coth(Pe) - \frac{1}{Pe} \quad (102) $$

where the Peclet number $Pe$ is given by

$$ Pe = \frac{E_1(\varphi) \cdot h}{2 \cdot E_2(\varphi)} \quad (103) $$

Integration of equation (99) in $]-\pi, \pi]$ leads to the ordinary differential equations system

$$ \underline{\Psi}^{*T} \underline{M} \dot{\underline{\Psi}} + \underline{\Psi}^{*T} \underline{G}_0 \underline{\Psi} + \underline{\Psi}^{*T} \underline{G}_1 \underline{\Psi} + \underline{\Psi}^{*T} \underline{G}_2 \underline{\Psi} = 0 \quad (104) $$

Eq.(104) can be rewritten introducing the matrix $\underline{G}$, with $\underline{G} = \underline{G}_0 + \underline{G}_1 + \underline{G}_2$, as

$$ \underline{M} \dot{\underline{\Psi}} + \underline{G} \underline{\Psi} = 0 \quad (105) $$

The system defined by the last equation under the normality condition, can be solved in an usual manner. The model reduction technique described in the previous sections can be applied for solving this ODE system using a very reduced number of degrees of freedom.
## 4.3 Results and Discussion

An homogenous shear test has been considered to illustrate the model reduction capabilities. An isotropic orientation distribution is assumed as initial orientation state. The velocity gradient is given by

$$ \mathrm{Grad}\underline{v} = \begin{pmatrix} 0 & \dot{\gamma} \\ 0 & 0 \end{pmatrix} \qquad (106) $$

with $\dot{\gamma} = 1$. Orientation space $\varphi \in [-\pi, \pi]$ has been divided in 72 finite elements. Basis enrichment operates when the norm of the residual related to the reduced model solution exceeds 0.001. Each enrichment is performed using the two first Krylov-subspaces. Moreover, a Karhunen-Loève decomposition allows to select the most significant basis functions. We consider all the eigenvectors associated with eigenvalues $\lambda_i > 10^{-8}\lambda_1$ ($\lambda_1$ being the highest eigenvalue).

<figure><img src="image_15.png" /><figcaption>Figure 14. Time-evolution of the distribution function in the start-up shearing test: $k = 0.9$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

The first case treated deal with a start-up shearing test during 10s where the steady state is almost reached. The factor $k$ is set to 0.9 and the diffusion coefficient $D_r$ is set to 0.1. Figure 14 shows the time evolution of the computed solution using the reduced model that consists of the 8 functions depicted in Figure 15. The evolution of the coefficients related to each function involved in the reduced model is depicted in Figure 16. The same problem has been solved using the finite element method involving 72 degrees of freedom. The error between the finite element and the reduced model solutions is shown in Figure 17. It can be clearly noticed the high accuracy provided by the reduced model solution despite the reduced number of degrees of freedom involved (8 instead of 72).

A second test is conducted with the same material and simulation parameters with the exception of the shear rate which is now time dependant. The shear rate is prescribed according to:

$$ \left\{ \begin{array}{l} \dot{\gamma} = +1 \text{ for time in } [0, 5[ \\ \dot{\gamma} = -1 \text{ for time in } [5, 10[ \\ \dot{\gamma} = +1 \text{ for time in } [10, 15[ \\ \dot{\gamma} = -1 \text{ for time in } [15, 20] \end{array} \right. \qquad (107) $$

Figure 18 shows the reduced model solutions at times $t=0, 0.5, 1, 2, 3, 5, 5.5, 6, 7, 8, 10, 10.5, 11, 12, 13, 15, 15.5, 16, 17, 18$ and $20$. The size of the reduced basis was found equal
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Reduced basis functions used for computing the solution depicted in Figure 14: $k = 0.9$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 16.** Evolution of the reduced model coefficients in the shear test: $k = 0.9$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

to 11 at the end of simulation (without any basis reduction). The approximation basis was reduced using the Karhunen-Loève decomposition for obtaining the 4 functions depicted in figure 19. The evolution of the coefficients related to the reduced approximation basis functions is shown in the Figure 20. For this calculation the basis enrichment has been performed only in the time interval [0, 10]. For all times greater than $t = 10$, calculation has been done without basis enrichment. Thus we can conclude that once reduced basis was built with a specific solicitation, the basis functions can reproduce accurately the solution for the same kind of solicitations.

In the third and last example the same kind of cyclic solicitation has been applied by changing the fiber parameter $k$ which is now set to 0.5. Figure 21 shows the reduced model solutions at different times, computed by using the approximation functions depicted in Figure 22. Finally, in Figure 23 the evolution of the different reduced basis coefficients is
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Comparison between the finite element and the reduced model solutions in the shear test: *k* = 0.9, *D*<sub>r</sub> = 0.1, ȳ = 1</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 18.** Time-evolution of the distribution function: *k* = 0.9, *D*<sub>r</sub> = 0.1, ȳ = +1/-1/+1/-1</figcaption></figure>

shown. In this case only 5 approximation functions are required in the simulation due to
the regularity of the searched solution related to *k* = 0.5. However, the FEM needs the
same number of shape functions in both cases. Thus, we can conclude that the number of
approximation functions involved in the reduced modelling depends only on the regularity
of the problem solution.

**5 CONCLUSIONS**

In this work we describe an algorithm which builds basis functions for function decomposi-
tion, thanks to an incremental method. This method involves a Karhunen-Loève expansion
of reduced state variables. This leads to an optimal number of basis functions. The pro-
posed algorithm is based on the APHR method. This is an adaptive strategy to build
<figure><img src="image_20.png" /><figcaption>**Figure 19.** The four most significant functions of the reduced basis involved in the reduced model: $k = 0.9, D_r = 0.1, \dot{\gamma} = +1/-1/+1/-1$</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Evolution of the reduced model coefficients: $k = 0.9, D_r = 0.1, \dot{\gamma} = +1/-1/+1/-1$</figcaption></figure>

reduced order model when the state evolution is implicitly defined by linear or non-linear governing equations. In case of known state evolutions the APHR method is an incremental Karhunen-Loève decomposition. A detailed formal example is given to illustrate a reduction of the number of basis functions. A reduced integration domain is proposed to illustrate the capability of the hyperreduction (reduction of the number of basis functions and reduction of the integration domain).

In the second part of this work we have analyzed the application of model reduction in transient problems defined in domains evolving in time. For this purpose we have solved some examples of uncoupled/coupled, linear/non-linear advection-diffusion equations, involving or not nodal redistributions (remeshing in the finite element context). We have
<figure><img src="image_22.png" /><figcaption>**Figure 21.** Time-evolution of the distribution function in the start-up shearing test: $k = 0.5$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 22.** Reduced model functions involved in the start-up shearing test: $k = 0.5$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

proved that the use of an “a priori” model reduction strategy in combination with a nodal updating according to the advection field and consequently with the domain evolution, avoids the necessity of field projections and stabilization for accounting advection effects. Obviously, a nodal updating without nodal redistribution requirements is only possible when one proceeds in the framework of meshless techniques, which are more flexible than the finite element method with respect to the relative position of nodes. Among the vast family of meshless techniques, we propose the use of the natural element method which, as previously described, is the only where the imposition of essential boundary conditions is as simple as in finite elements.

Finally a kinetic theory model has been solved to illustrate the capabilities of the proposed technique for simulating complex flows involving fluids with microstructure (polymer solutions, polymer melts or suspensions, among many others).
<figure><img src="image_24.png" /><figcaption>Figure 23. Evolution of the reduced model coefficients: $k = 0.5$, $D_r = 0.1$, $\dot{\gamma} = 1$</figcaption></figure>

REFERENCES

1 A. Ammar and F. Chinesta (2004). 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 Engng.* Springer, **43**, 1-16.

2 G. P. Brooks and J. M. Powers (2002). A Karhunen-Loève Galerkin Technique with Shock Fitting for Optimization of a Blunt Body Geometry. In *38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and exhibit. Indianapolis, IN*.

3 F. Chinesta, G. Chaidron and A. Poitou (2003). 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), 97-125.

4 E. Cueto, M. Doblaré and L. Gracia (2000). Imposing essential boundary conditions in the Natural Element Method by means of density-scaled $\alpha$-shapes. *International Journal for Numerical Methods in Engineering*, **49**(4), 519-546.

5 E. Cueto, N. Sukumar, B. Calvo, J. Cegoñino and M. Doblaré (2003). Overview and recent advances in Natural Neighbour Galerkin methods. *Archives of Computational Methods in Engineering*, **10**(4), 307-384.

6 M. Fahl (2001). Computation of POD Basis Functions for Fluid Flows with Lanczos Methods. *Mathematical and Computer Modeling*, **34**, 91-107.

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

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

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

10 M. M. Loève (1963). *Probability theory*. The University Series in Higher Mathematics, 3rd Ed. Van Nostrand, Princeton, NJ.
11 E. N. Lorenz (1956). *Empirical Orthogonal Functions and Statistical Weather Prediction*. MIT, Departement of Meteorology, Scientific Report N1, Statistical Forecasting Project.

12 M. A. Martínez, E. Cueto, M. Doblaré and F. Chinesta (2003). Natural Element meshless simulation of injection processes involving short fiber suspensions. *Journal of Non-Newtonian Fluid Mechanics*, **115**, 51–78.

13 H. M. Park and D.H. Cho (1996). The Use of the Karhunen-Loève Decomposition for the Modelling of Distributed Parameter Systems. *Chem. Engineer. Science*, **51**, 81–98.

14 D. Ryckelynck (2005). A priori hypereduction method: an adaptive approach. *Journal of Computational Physics*, **202**, 346–366.

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

16 N. Sukumar, B. Moran and T. Belytschko (1998). The Natural Element Method in Solid Mechanics. *International Journal for Numerical Methods in Engineering*, **43**(5), 839–887.