Generalization of the Neville-Aitken Interpolation Algorithm on Grassmann Manifolds : Applications to Reduced Order Model

TL;DR

This paper extends Neville-Aitken interpolation to Grassmann manifolds for parametric PDEs, enabling efficient reduced order modeling with demonstrated accuracy and computational benefits across multiple CFD applications.

Contribution

It introduces a novel recursive interpolation method on Grassmann manifolds using geodesic barycenters, applicable to reduced order models in parametric PDEs.

Findings

Accurate interpolation of solution subspaces in CFD problems.

Reduced computational time for parametric simulations.

Effective application across diverse fluid dynamics scenarios.

Abstract

The interpolation on Grassmann manifolds in the framework of parametric evolution partial differential equations is presented. Interpolation points on the Grassmann manifold are the subspaces spanned by the POD bases of the available solutions corresponding to the chosen parameter values. The well-known Neville-Aitken's algorithm is extended to Grassmann manifold, where interpolation is performed in a recursive way via the geodesic barycenter of two points. The performances of the proposed method are illustrated through three independent CFD applications, namely: the Von Karman vortex shedding street, the lid-driven cavity with inflow and the flow induced by a rotating solid. The obtained numerical simulations are pertinent both in terms of the accuracy of results and the time computation.

Tables

Table 1. Table 1 : Description of the cases for the application of section  4.1 (Von Karman vortex shedding street). The set of sampling parameters for which the simulation of the full order model is performed is Λ Λ \Lambda . The interpolation is performed at the new parameter λ ⋆ superscript 𝜆 ⋆ \lambda^{\star} .

Case	Sampling $\Lambda$	Target ${\lambda }^{\star }$
$1$	$(100,120,130,160,170,200)$	$110$
$2$	$(100,160,170,180,200)$	$110$
$3$	$(100,120,130,140,200)$	$190$

Table 2. Table 2 : Results for the projection error defined in ( 17 ) for the application of section  4.1 (Von Karman vortex shedding street). Reference refers to the basis built directly from the snapshots obtained by the simulation of the HDM at the target parameter λ ⋆ superscript 𝜆 ⋆ \lambda^{\star} (no interpolation). Neville , Amsallem and Standard refer to the bases obtained from the interpolation of the sampling sets defined in table  1 (no simulation of the HDM).

Method	Case 1	Case 2	Case 3
Reference	6.89e-05	6.89e-05	3.10e-04
Neville	1.03e-03	7.95e-03	1.64e-03
Amsallem	2.54e-03	9.12e-02	6.68e-03
Standard	8.51e-01	2.13e-01	1.41e-01

Table 3. Table 3 : Results for the error associated with the dynamical systems as defined in ( 18 ) for the application of section  4.1 (Von Karman vortex shedding street).

Method	Case 1	Case 2	Case 3
Reference	4.41e-04	4.41e-04	1.37e-03
Neville	1.74e-03	1.35e-02	4.46e-03
Amsallem	3.42e-03	1.05e-01	1.31e-02
Standard	9.36e-01	3.06e-01	1.60e-01

Table 4. Table 4 : Description of the cases for the application of section  4.2 (lid-driven cavity with inflow). The parameter is the ratio of the Reynolds number associated with the lid velocity and that associated with the inflow velocity as defined in ( 20 ). The set of sampling parameters for which the simulation of the full order model is performed is Λ Λ \Lambda . The interpolation is performed for the new parameter  λ ⋆ superscript 𝜆 ⋆ \lambda^{\star} .

Case	Sampling $\Lambda$	Target ${\lambda }^{\star }$
1	(0.0, 0.2, 0.4, 0.6, 0.8, 1.0)	0.5
2	(0.0, 0.2, 0.3, 0.4, 0.9, 1.0)	0.5
3	(0.0, 0.2, 0.3, 0.4, 0.9, 1.0)	0.7
4	(0.0, 0.3, 0.4, 0.9, 1.0)	0.1

Table 5. Table 5 : Results for the projection error as defined in ( 17 ) for the application of section  4.2 (lid-driven cavity with inflow). Reference refers to the basis built directly from the snapshots obtained by the simulation of the HDM at the target parameter λ ⋆ superscript 𝜆 ⋆ \lambda^{\star} (no interpolation). Neville and Amsallem refer to the bases obtained from the interpolation of the sampling sets defined in table  4 (no simulation of the HDM).

Method	Case 1	Case 2	Case 3	Case 4
Reference	2.98e-06	2.98e-06	4.71e-05	1.90e-06
Neville	2.7e-02	2.2e-01	2.0e-01	1.6e-01
Amsallem	2.8e-02	2.1e-01	4.6e-01	4.8e-01

Table 6. Table 6 : Results for the error associated with the dynamical systems as defined in ( 18 ) for the application of section  4.2 (lid-driven cavity with inflow).

Method	Case 1	Case 2	Case 3	Case 4
Reference	1.56e-03	1.56e-03	4.54e-03	5.13e-04
Neville	4.6e-02	3.5e-01	5.1e-01	2.7e-01
Amsallem	5.1e-02	3.5e-01	9.7e-01	6.2e-01

Table 7. Table 7 : Description of the cases for the application of section  4.3 (flow induced by a rotating body). The parameter is the Reynolds number as defined in ( 21 ). The set of sampling parameters for which the simulation of the full order model is performed is Λ Λ \Lambda . The interpolation is performed for the new parameter  λ ⋆ superscript 𝜆 ⋆ \lambda^{\star} .

Case	Sampling $\Lambda$	Target ${\lambda }^{\star }$
$1$	$(500,1000,1500,2500,3000)$	$2800$
$2$	$(500,1000,1500,1800,2500,2800,3000)$	$800$
$3$	$(500,1500,3000)$	$2800$

Table 8. Table 8 : Results for the projection error as defined in ( 17 ) for the application of section  4.3 (flow induced by a rotating body).

Method	Case 1	Case 2	Case 3
Reference	6.493e-04	2.747e-04	6.493e-04
Neville	1.874e-03	2.457e-03	2.141e-03
Amsallem	1.880e-03	2.457e-03	2.143e-03

Table 9. Table 9 : Results for the error associated with the dynamical systems as defined in ( 18 ) for the application of section  4.3 (flow induced by a rotating body).

Method	Case 1	Case 2	Case 3
Reference	2.964e-02	3.029e-03	2.964e-02
Neville	2.108e-02	1.817e-02	1.474e-02
Amsallem	2.189e-02	1.818e-02	1.475e-02

Equations

$$X\sim Y,\text{if and only if}\exists A\in G{L}_m(\mathbb{R})\text{such that :}X=YA,$$

$$G_m({\mathbb{R}}^n)=\{\overline{X}:X\in {\mathbb{R}}_{\ast }^{n\times m}\},$$

$$\overline{X}=\{XA:A\in G{L}_m(\mathbb{R})\}.$$

$$\pi :{\mathbb{R}}_{\ast }^{n\times m}\to G_m({\mathbb{R}}^n),$$

$$g_{{}_{XA}}(Y_{1}A,Y_{2}A)=g_{{}_X}(Y_1,Y_2),$$

$$g_{{}_X}(Y_1,Y_2)=tr[{(}^{t}XX{)}^{-1}{}^t{Y}_1{Y}_2],$$

$$d_X\pi :H_{{}_X}⟶T_{\overline{X}}{G}_m({\mathbb{R}}^n)$$

$$(d_X\pi {)}^{-1}(v){\left({}^{t}XX\right)}^{1/2}=U\Sigma {}^{t}V,$$

$$\gamma (s)=\overline{X{(}^{t}XX{)}^{-1/2}V\cos (s\Sigma )+U\sin (s\Sigma )}.$$

$$X{(}^{t}XX{)}^{-1/2}V\cos (\Sigma )+U\sin (\Sigma ).$$

$${\gamma }_{{}_{\bar{X},\bar{Y}}}(s)=\overline{X\tilde{V}\cos \left(s{\tan }^{-1}\left(\tilde{\Sigma }\right)\right)+\tilde{U}\sin \left(s{\tan }^{-1}\left(\tilde{\Sigma }\right)\right)},$$

$$d(\overline{X},\overline{Y})=\Vert {\tan }^{-1}\left(\tilde{\Sigma }\right){\Vert }_{{}_F},$$

$$r_{{}_I}\left(\overline{X}\right)=\sup \left\{\rho >0:{\exp }_{\overline{X}}:B(0_{\overline{X}},\rho )\to {\exp }_{\overline{X}}(B(0_{\overline{X}},\rho ))\text{is a diffeomorphism}\right\},$$

$$r_{{}_I}=\inf \{r_{{}_I}\left(\overline{X}\right):\overline{X}\in G_m({\mathbb{R}}^n)\}.$$

$$U_{\alpha }{\tan }^{-1}({\Sigma }_{\alpha }){}^t{V}_{\alpha },$$

$$U_{\alpha }{\Sigma }_{\alpha }^t{V}_{\alpha }=[Y{(}^{t}XY{)}^{-1}{(}^{t}XX)-X]{(}^{t}XX{)}^{-1/2}.$$

$$L[{\lambda }_0,\cdots ,{\lambda }_N;y_0,\cdots ,y_N]({\lambda }_k)=y_k,\text{for every}k\in \{0,1,\cdots ,N\}.$$

$${\alpha }_{{}_{[i,j]}}(\lambda )=\frac{\lambda -{\lambda }_i}{{\lambda }_j-{\lambda }_i},\forall i,j\in \{0,1,\cdots \}.$$

$$L[{\lambda }_0,\cdots ,{\lambda }_N;y_0,\cdots ,y_N](\lambda )=$$

$$+{\alpha }_{{}_{[0,N]}}(\lambda )$$

$$L[{\lambda }_k;y_k](\lambda )=y_k\text{for every}k\in \{0,1,\cdots ,N\}.$$

$$L[{\lambda }_0,{\lambda }_1;y_0,y_1](\lambda )$$

$$L[{\lambda }_1,{\lambda }_2;y_1,y_2](\lambda )$$

$$L[{\lambda }_0,{\lambda }_1,{\lambda }_2;y_0,y_1,y_2](\lambda )$$

$$\text{(H)}\exists y\in G_m({\mathbb{R}}^n):\forall k\in \{0,1,\cdots ,N\},y_k\in B_d(y,\pi /2).$$

$$\Gamma [{\lambda }_i,{\lambda }_j;y_i,y_j]:={\gamma }_{{}_{[y_i,y_j]}}\circ {\alpha }_{{}_{[i,j]}}$$

$$\Gamma [{\lambda }_i,{\lambda }_j;y_i,y_j](\lambda )=\underset{y\in G_m({\mathbb{R}}^n)}{\mathrm{argmin}}J(y),$$

$$J(y)=\frac{1}{2}\left[\left(1-{\alpha }_{{}_{[i,j]}}(\lambda )\right)d^2(y,y_i)+{\alpha }_{{}_{[i,j]}}(\lambda )d^2(y,y_j)\right],$$

$$\left(1-{\alpha }_{{}_{[i,j]}}(\lambda )\right){\exp }_{{}_{\Gamma [{\lambda }_i,{\lambda }_j;y_i,y_j](\lambda )}}^{-1}(y_i)+{\alpha }_{{}_{[i,j]}}(\lambda ){\exp }_{{}_{\Gamma [{\lambda }_i,{\lambda }_j;y_i,y_j](\lambda )}}^{-1}(y_j)=0.$$

$${\gamma }_1(t)$$

$${\gamma }_2(t)$$

$${\gamma }_1^{\prime }(0)$$

Full text

Keywords : Neville interpolation algorithm, Grassmann Manifold, Lagrange interpolation

Mathematics Subject Classification (2010) : 65K05, 65D18, 65F15.

1 Introduction

In the present manuscript, we are interested with interpolation on Grassmann manifolds in the framework of parametric evolution partial differential equations. A relevant example in this topic is the Navier-Stokes equations where the varying parameter $\lambda$ is the Reynolds number. We confine ourselves in the situation where the solution of the parametric partial differential equation is computed for a set $\{\lambda_{0},\,\lambda_{1},\cdots,\lambda_{N}\}$ of parameter values and the Proper Orthogonal Decomposition (POD) is used to construct the associated bases $\Phi(\lambda_{0}),\,\Phi(\lambda_{1}),\cdots\Phi(\lambda_{N})$. These solutions are in general computed with high computational costs discretization type methods as Finite Elements, Finite Volumes, Discontinuous Galerkin, etc.

For a new value of the parameter $\lambda$, instead of recomputing the solution, we can use the already computed POD bases to carry out an interpolation. A natural question is how to interpolate a set of given bases. The purpose of the present work is to develop low-dimensional parametric methods to compute an approximation of the solution, by interpolation of POD bases, for a new parameter value.

In various problems in fluid mechanics or fluid structure interaction, the POD bases are computed for a given Reynolds number or for other flow parameters. Theses bases are then used to build a reduced model to predict the flow for other values of these parameters. The parametric domain of validity of such reduced models has to be described. In the framework of quasilinear parabolic problems, this parametric domain of validity was studied according to the solution’s regularity and the number of modes retained to build the reduced model by Akkari et al. [5]. These results were then extended to Navier-Stokes equations [3] and Burgers equation [4]. Numerical results on parametric sensitivity for fluid structure interaction problems were carried out in [16]. In [17], the authors introduced non-intrusive strategies to provide vademecums intended to real-time computations for nonlinear thermo-mechanical problems.

In control theory, the numerical costs (CPU and memory) associated with the adjoint equation-based methods used to solve the underlying optimization problems are so important that the three-dimensional Navier-Stokes equations are rarely studied.

In [14], the POD is applied to solve open-loop and closed-loop optimal control problems for the Burgers equation. Similar methods are used to control laser surface hardening, where the state equations are a semilinear heat equation coupled with an ordinary differential equation [21]. Phase-field control problems are also investigated by [29]. In medical imaging [9, 26], parametric POD bases are studied when time-interval variation, spatial domain variation or parameter variation may occur.

One classical way to perform such approximations is the use of Reduced Basis (RB) methods. RB-methods consist in the approximation of the manifold solution $\mathcal{M}\subset H$, i.e., the set of parametric solutions, by a low-dimensional subspace $H_{m}\subset H$, where $H$ is the Hilbert space in which the functional framework of the problem is well-posed. In this topic, one popular way is the construction of such a subspace via the snapshots $\left(u(\lambda_{k})\right)_{0\leq k\leq{N}}$ as, for example, $H_{m}={\rm span}\{u(\lambda_{k})\;:\;0\leq k\leq{N}\}$ and use Galerkin type techniques to compute an approximation $u_{{}_{m}}(\lambda)\in H_{m}$ of $u(\lambda)\in\mathcal{M}$. We refer the interested reader to [25, 19] and the references therein.

Another way to approximate the parametric solution $u$ for a new parameter value $\lambda\in\mathbb{R}^{P}$ from given data $\left(u(\lambda_{k})\right)_{0\leq k\leq{N}}$ via POD techniques is the interpolation on Grassmann manifolds. We describe briefly this method and refer the reader to the work of Amsallem and Farhat [6] for more details. This procedure can be summarized as follows: for each $k\in\{0,1,\cdots,{N}\}$, consider the rank-$m$ POD basis $\left(\Phi_{{}_{i}}(\lambda_{k})\right)_{1\leq i\leq m}$ and the corresponding spanned subspace $H_{m}^{(k)}={\rm Span}\left\{\Phi_{{}_{i}}(\lambda_{k})\;:\;1\leq i\leq m\right\}\subset H$, (in [6], the space $H$ has a large but finite dimension). Now, denote by $G_{m}(H)$ the Grassmann manifold consisting of all $m$-dimensional subspaces of $H$. Then the $m$ subspaces $H_{m}^{(k)}$ can be seen as points on this manifold $G_{m}(H)$. For a new parameter value $\lambda$, the question is how to compute a good approximation of the rank-${m}$ POD basis associated with the solution $u(\lambda)$, or in a more relevant way, how to find its corresponding point in the Grassmann manifold $G_{m}(H)$? Indeed, it is known from the work of Amsallem and Farhat [6] that the appropriate objects to interpolate are not the POD bases but the underlying spanned subspaces. To answer this question, the idea developed in [6] is to choose a reference parameter value $\lambda_{r}\in\{\lambda_{k}\;:\;0\leq k\leq{N}\}$, to consider the POD subspace $R:=H_{m}^{(r)}$ as a reference point on the Grassmann manifold $G_{m}(H)$. Let us denote by $T_{{}_{R}}G_{m}(H)$ the tangent space to the manifold $G_{m}(H)$ at the point $H_{m}^{(r)}$. Then, the logarithm mapping at the reference point $H_{m}^{(r)}$ denoted by $\log^{(r)}\;:\;G_{m}(H)\longrightarrow T_{{}_{R}}G_{m}(H)$ allows to map (in a local way) all points $\left(H_{m}^{(k)}\right)_{0\leq k\leq N}$ on the tangent space $T_{{}_{R}}G_{m}(H)$ where one can perform an interpolation (see Section 2 for more details). Finally, one returns to the manifold $G_{m}(H)$ via the (reciprocal) exponential mapping from the tangent space $T_{{}_{R}}G_{m}(H)$ to get the desired approximation of $H_{m}^{(\lambda)}$ on $G_{m}(H)$.

The disadvantage of the method developed in [6] is its dependence on the reference point and the lack of a strategy for choosing this reference point.

In a recent work [18], the authors developed the Grassmann Inverse Distance Weighting (G-IDW) as an extension of the well-known Inverse Distance Weighting (IDW) method to Grassmann manifolds. It consists in the minimization of a quadratic function in the geodesic distance to interpolation points, with appropriate weights. This G-IDW method does not require a reference point and yields relevant results. As a counterpart, it is iterative.

In the present work, we propose the extension of the Neville-Aitken’s algorithm which computes the Lagrange interpolation polynomial in a recursive way from the interpolation of two points. Replacing straight lines by geodesics, we extend this algorithm to Grassmann manifolds. The obtained numerical simulations are excellent both in terms of accuracy of the results and the computation time.

The present paper is organized as follows: In Section 2, the Grassmann manifold and the underlying tools used in the paper are briefly presented. Section 3 is dedicated to the introduction of the Neville-Aitken’s algorithm in the framework of vector spaces and its extension to Grassmann manifolds. Finally, the performances of the proposed method are illustrated through three independent CFD applications, namely: the Von Karman vortex shedding street, the lid-driven cavity with inflow and the flow induced by a rotating solid. Our results are compared to those obtained by the exact POD method and to the interpolation approach developed by Amsallem et al. [6].

2 The Grassmann manifold

Let $m$, $n$ be two positive integers with $m\leq n$. The set of all $\mathbb{R}$-subspaces of $\mathbb{R}^{n}$ of dimension $m$ is a differentiable manifold of dimension $m\times(n-m)$, called the Grassmann manifold and usually denoted by $G_{m}(\mathbb{R}^{n})$. Let $\mathbb{R}^{n\times m}_{*}$ be the set of all matrices of size $n\times m$, whose column vectors are linearly independent; one defines on this set the following equivalence relation

[TABLE]

where $GL_{m}(\mathbb{R})$ denotes the linear group of degree $m$, i.e., the set of $m\times m$ invertible matrices. The Grassmann manifold can be realized as the set of such equivalence classes. More precisely,

[TABLE]

where

[TABLE]

The Grassmann manifold is a compact topological space and, thanks to the canonical projection denoted by

[TABLE]

it can be endowed by a differentiable manifold structure so that $\pi$ is a submersion ([23]).

2.1 Exponential and logarithm mappings, injectivity radius

On each point $X\in\mathbb{R}^{n\times m}_{*}$, one constructs a metric $g_{{}_{X}}$ which is right-invariant by $GL_{m}(\mathbb{R})$, i.e.,

[TABLE]

defined by

[TABLE]

for every $Y_{1},Y_{2}\in\mathbb{R}^{n\times m}$ and $A\in GL_{m}(\mathbb{R})$. This metric induces a Riemannian metric on $G_{m}(\mathbb{R}^{n})$ through the submersion $\pi$. Moreover, one constructs a Riemannian connection structure on $G_{m}(\mathbb{R}^{n})$, which allows the definition of geodesics and the geodesic *exponential * mapping on $G_{m}(\mathbb{R}^{n})$ [1].

In what follows, the distance on $G_{m}(\mathbb{R}^{n})$, induced by the metric $g$, will be denoted by $d$, and the ball centered at $y$ of radius $r>0$ in $G_{m}(\mathbb{R}^{n})$ will be denoted by $B_{{}_{d}}(y,r)$.

For every $X\in\mathbb{R}^{n\times m}_{*}$, we introduce $H_{{}_{X}}=\left[T_{{}_{X}}\pi^{-1}\left(\overline{X}\right)\right]^{\perp}$, the orthogonal space to the tangent space at $X$ of the fiber $\pi^{-1}\left(\overline{X}\right)$, where the orthogonality is with respect to scalar product induced by the metric $g$. Then the mapping

[TABLE]

is an isomorphism. Now, given a vector $v\in T_{\overline{X}}\,G_{m}(\mathbb{R}^{n})$, if we perform the singular value decomposition (SVD) of the matrix $(d_{X}\pi)^{-1}(v)\,\left({}^{t}XX\right)^{1/2}$ as

[TABLE]

the geodesic of initial point $\gamma(0)=\overline{X}$ and initial velocity $\gamma^{\prime}(0)=v$ can be written in the form

[TABLE]

The geodesic exponential mapping is defined by $\exp_{\overline{X}}(v)=\gamma(1)$ and has a matrix representation given by:

[TABLE]

On the other hand, since $d_{0_{{}_{\overline{X}}}}\exp_{\overline{X}}=I$, then the inverse function theorem implies that $\exp_{\overline{X}}$ is a local diffeomorphism at $0_{{}_{\overline{X}}}$, the null vector of $T_{\overline{X}}G_{m}\left(\mathbb{R}^{n}\right)$. In the same way, we can define the geodesic between two points $\overline{X},\overline{Y}$ in $G_{m}(\mathbb{R}^{n})$, with ${}^{t}XY\in GL_{m}(\mathbb{R})$, by

[TABLE]

where $\widetilde{U}\>\widetilde{\Sigma}\>^{t}\widetilde{V}=Y(^{t}X\>Y)^{-1}-X$ is a SVD.

Notice that

[TABLE]

where $\|\,\cdot\,\|_{{}_{F}}$ denotes the Frobenius norm on the space of matrices.

Definition 1** (injectivity radius).**

Let $\overline{X}\in G_{m}(\mathbb{R}^{n})$. The injectivity radius of $(G_{m}(\mathbb{R}^{n}),g)$ at $\overline{X}$ is defined by

[TABLE]

and the injectivity radius of $G_{m}(\mathbb{R}^{n})$ is defined by

[TABLE]

The value of the injectivity radius of $G_{m}(\mathbb{R}^{n})$ is given by the following result.

Theorem 1**.**

[31]** Let $n,m$ be two positive integers such that $min\{m,n-m\}\geq 2$. Then the injectivity radius of $G_{m}(\mathbb{R}^{n})$ is given by $r_{{}_{I}}=\pi/2$.

The proof of Theorem 1 can be found in [13]. This Theorem gives the radius of the ball where the reciprocal mapping of the geodesic exponential $\exp_{\overline{X}}^{-1}$ (or $\log_{\overline{X}}$) is defined. More precisely, for $\overline{Y}\in G_{m}(\mathbb{R}^{n})$, with $\displaystyle d(\overline{X},\overline{Y})<\frac{\pi}{2}$, then $\exp_{\overline{X}}^{-1}(\overline{Y})$ is represented by the matrix

[TABLE]

with the singular value decomposition

[TABLE]

In the next section, we will introduce an extension of the well-known Neville-Aitken algorithm [2, 24] to the Grassmann manifold.

3 Neville-Aitken’s algorithm for Lagrange interpolation

In this section, we first recall the classical Neville-Aitken’s algorithm for vector-valued functions, and second introduce our extension to the Grassmann manifolds.

3.1 Recall of the classical Neville-Aitken’s algorithm

Consider a function $f\,:\,\Lambda\subset\mathbb{R}\longrightarrow V$, where $V$ is a vector space, and a subset $\{\lambda_{0},\lambda_{1},\cdots\}\subset\Lambda$, with $f(\lambda_{k})=y_{k}\in V$, for every $\lambda_{k}$.

Given an integer $N\geq 1$, the Lagrange interpolation polynomial of $f$ at $\{\lambda_{0},\lambda_{1},\cdots,\lambda_{N}\}$ will be denoted by $L[\lambda_{0},\cdots,\lambda_{N};y_{0},\cdots,y_{N}]$. It is the unique polynomial of degree less or equal to $N$, satisfying:

[TABLE]

Let us consider the affine functions $\alpha_{{}_{[i,j]}}$ defined on $\Lambda$ by

[TABLE]

Then, the Neville-Aitken’s algorithm can be summarized by the following induction formula:

[TABLE]

with the initialization

[TABLE]

That is, for a fixed $\lambda\in\Lambda$, the quantity $L[\lambda_{0},\cdots,\lambda_{N};y_{0},\cdots,y_{N}](\lambda)$ is the barycenter of

$L[\lambda_{0},\cdots,\lambda_{N-1};y_{0},\cdots,y_{N-1}](\lambda)$ and $L[\lambda_{1},\cdots,\lambda_{N};y_{1},\cdots,y_{N}](\lambda)$, weighted by $(1-\alpha_{{}_{[0,N]}}(\lambda))$ and $\alpha_{{}_{[0,N]}}(\lambda)$, respectively. Therefore, Neville-Aitken’s algorithm reduces the Lagrange interpolation into recursive barycenter computations of two points at each step, with appropriate weights. Thus, for a fixed value of the parameter $\lambda$, we get

[TABLE]

with ${\rm deg}(L[\lambda_{0},\lambda_{1};y_{0},y_{1}])\leq 1,$ ${\rm deg}(L[\lambda_{1},\lambda_{2};y_{1},y_{2}])\leq 1,$ ${\rm deg}(L[\lambda_{0},\lambda_{1},\lambda_{2};y_{0},y_{1},y_{2}])\leq 2,$ and so on to obtain the Lagrange interpolation polynomial $L[\lambda_{0},\cdots,\lambda_{N};y_{0},\cdots,y_{N}]$, whose degree is less or equal to $N$.

Remark 1**.**

For vector-valued functions, the Lagrange interpolation polynomial is invariant by permutations of interpolation points. That is, if $\sigma$ is a permutation of the set $\{\lambda_{0},\lambda_{1},\cdots,\lambda_{N}\}$, then $L[\lambda_{\sigma{(0)}},\cdots,\lambda_{\sigma{(N)}};y_{\sigma{(0)}},\cdots,y_{\sigma{(N)}}]=L[\lambda_{0},\cdots,\lambda_{N};y_{0},\cdots,y_{N}]$.

3.2 Extension of the Neville-Aitken’s algorithm to Grassmann manifolds

As pointed before, if we are able to define the notion of barycenter of two points on the Grassmann manifold, then we will be able to carry out interpolation on it, via Neville-Aitken’s algorithm.

More precisely, consider a function $f\,:\,\Lambda\subset\mathbb{R}\longrightarrow G_{m}(\mathbb{R}^{n})$ such that $y_{i}:=f(\lambda_{i})$ are given, for $i\in\{0,1,\cdots\}$. We will describe one way to handle with the interpolation of $f$ at $\{\lambda_{0},\lambda_{1},\cdots,\lambda_{N}\}$, for any $N\geq 1$.

Remark 2**.**

Since our interpolation extension will be based on geodesics on the Grassmann manifold $G_{m}(\mathbb{R}^{n})$, we will assume hereafter that all the points $\left(y_{i}\right)_{0\leq i\leq N}$ belong to a certain ball of radius $\frac{\pi}{2}$ in $G_{m}(\mathbb{R}^{n})$. That is

[TABLE]

The assumption (H) will ensure the uniqueness of geodesics connecting each pair of interpolation points [31].

At first, consider two (distinct) values $\lambda_{i}$ and $\lambda_{j}$ whose images $y_{i}$ and $y_{j}$ belong to the Grassmann manifold $G_{m}(\mathbb{R}^{n})$. Let $\gamma_{{}_{[y_{i},y_{j}]}}\,:\,[0,1]\longrightarrow G_{m}(\mathbb{R}^{n})$ be the unique geodesic on $G_{m}(\mathbb{R}^{n})$ satisfying $\gamma_{{}_{[y_{i},y_{j}]}}(0)=y_{i}$ and $\gamma_{{}_{[y_{i},y_{j}]}}(1)=y_{j}$. Then, the map defined on $[\lambda_{i},\lambda_{j}]$ by

[TABLE]

is the unique geodesic on $G_{m}(\mathbb{R}^{n})$ satisfying $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda_{i})=y_{i}$ and $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda_{j})=y_{j}$, where $\alpha_{{}_{[i,j]}}$ is defined by (5). Indeed, geodesics are invariant by affine transformations and thanks to equation (2) a matrix representation of the geodesic connecting two points on the manifold is known. This matrix representation will be used in the algorithm 1.

Definition 2**.**

For a given $\lambda\in[\lambda_{i},\lambda_{j}]$, the point $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda)\in G_{m}(\mathbb{R}^{n})$ is called the barycenter (in the sense of geodesics) on $G_{m}(\mathbb{R}^{n})$ of   $y_{i}$ and $y_{j}$, weighted by $\left(1-\alpha_{{}_{[i,j]}}(\lambda)\right)$ and $\alpha_{{}_{[i,j]}}(\lambda)$ respectively.

The barycenter $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}]$ provides an interpolation on $G_{m}(\mathbb{R}^{n})$ of the function $f$ at $\lambda_{i}$ and $\lambda_{j}$. Moreover, it coincides with the barycenter in the sense of Karcher [12]. More precisely, we have the following result.

Proposition 1**.**

Let $\lambda_{i}<\lambda_{j}$ be two parameter values in $\{\lambda_{0},\lambda_{1},\cdots\}$ and $y_{i}$, $y_{j}$ be two points on the Grassmann manifold $G_{m}(\mathbb{R}^{n})$. Then

[TABLE]

where

[TABLE]

and $d$ denotes the geodesic distance on $G_{m}(\mathbb{R}^{n})$.

Proof. Applying [18, Theorem 4.1], it suffices to show that $\nabla J(\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda))=0$, or equivalently:

[TABLE]

To lighten the notations, we set $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda)=\widehat{y}$. We introduce the two geodesic paths on $\gamma_{1}$, and $\gamma_{2}$ defined from the interval $[0,1]$ to $G_{m}(\mathbb{R}^{n})$ by

[TABLE]

from $\widehat{y}$ to $y_{i}$ and from $\widehat{y}$ to $y_{j}$, respectively. Therefore, we get

[TABLE]

Whence, combining (11) and (12), we obtain

[TABLE]

which ends the proof. $\Box$

In the same way, we can extend the analogous of (6) to determine the interpolation $\Gamma[\lambda_{i},\lambda_{j},\lambda_{k}\,;\,y_{i},y_{j},y_{k}]$ of $f$ at $\{\lambda_{i},\lambda_{j},\lambda_{k}\}$. Indeed, we compute at first $\Gamma[\lambda_{i},\lambda_{j}\,;\,y_{i},y_{j}](\lambda)$ and $\Gamma[\lambda_{j},\lambda_{k}\,;\,y_{j},y_{k}](\lambda)$ by (9) and then

[TABLE]

and so on to obtain the interpolation $\Gamma[\lambda_{0},\cdots,\lambda_{N}\,;\,y_{0},\cdots,y_{N}]$ of $f$ at $\{\lambda_{0},\cdots,\lambda_{N}\}$, that is :

[TABLE]

Therefore, the resulting Neville-Aitken algorithm on Grassmann manifolds can be summarized as follows:

4 Applications

The performances of the proposed method are illustrated in this section through three CFD applications, namely (i) the Von Karman vortex shedding street, (ii) the lid-driven cavity with inflow and (iii) the flow induced by a rotating solid. For each application, a relevant dimensionless parameter $\lambda$ is selected and the test procedure is as follows.

We derive a High-Dimensional Model (HDM) by a standard finite-element discretization of the governing equations. 2. 2.

We build $N_{\lambda}$ POD bases by applying the method of snapshots [27] to the solution arising from the simulation of the HDM for a set of $N_{\lambda}$ sampling parameters $(\lambda_{i})_{i=0}^{N_{\lambda}}$. 3. 3.

We interpolate the sets of POD bases at a given new parameter $\lambda^{\star}$ by (i) the proposed method (labeled Neville) and (ii) that proposed in [6] (labeled Amsallem). In the latter case, the reference point of the method $\Phi(\lambda_{r})$ is chosen as the point in the sample bases closest to the target parameter: $\lambda_{r}=\lambda_{i}$ with $i=\underset{i\in\{0,\cdots,N_{\lambda}\}}{\operatorname*{argmin}}|\lambda^{\star}-\lambda_{i}|$. 4. 4.

We build a Reduced Order Model (POD-ROM) by a Gelerkin projection of the governing equations on the POD bases obtained from both interpolation methods. 5. 5.

The results are compared with the POD-ROM built from the reference POD basis computed directly from the snapshots of the HDM for $\lambda^{\star}$ (labeled Reference).

Details on the POD can be found in the reference [20] and in [30, 8]. All the numerical tests have been performed using the Python/C++ finite element library DOLFIN [15] on a computer equipped with a processor5552 sockets, 8 cores for each socket, 2 threads for each core, cadenced at 2.10GHz with a cache of 20MB. Intel Xeon E5-2620 v4 and 64Go of RAM.

4.1 Von Karman vortex shedding street

Here, we consider the planar flow of an incompressible newtonian fluid around a solid disk in a channel. The configuration is shown in figure 1. The HDM is derived by a standard variational formulation of the adimensional Navier-Stokes equations over the Taylor-Hood finite-elements space ($\mathrm{P2}$ for velocity and $\mathrm{P1}$ for the pressure). Here, the dimensionless parameter is the Reynolds number $\lambda={\frac{\rho\,v_{\infty}\,\ell}{\eta}}$ where $\rho=1$kg.m*-3* is the fluid density, $v_{\infty}=1$m.s*-1* is the magnitude of the velocity at inflow, $\ell=1$m is the disk diameter and $\eta$ is the dynamic viscosity. The test cases (i.e. the set of sampling parameters and the target parameter) are given in table 1.

For each sampling parameter $\lambda\in\Lambda$, the simulation is performed by a standard predictor/corrector scheme (see e.g. [11, §7.3]) and a set of ${N_{T}=200}$ regularly spaced snapshots of both the velocity \big{(}\mathbf{v}(\mathbf{x},t_{i},\lambda)\big{)}_{i=1}^{N_{T}} and the pressure \big{(}p(\mathbf{x},t_{i},\lambda)\big{)}_{i=1}^{N_{T}} are selected over a period $T={20}$s which does not include the transient period. To ensure homogeneous boundary conditions for the POD modes, both the velocity and the pressure are split into a mean field $\overline{\mathbf{v}}(\mathbf{x},\lambda)=\frac{1}{N_{T}}\sum_{i=1}^{N_{T}}\mathbf{v}(\mathbf{x},t_{i},\lambda)$ and a fluctuating field:

[TABLE]

Then, we compute the POD basis \Phi(\mathbf{x},\lambda)=\big{(}\phi_{m}(\mathbf{x},\lambda)\big{)}_{m=1}^{N_{T}} (respectively \Psi(\mathbf{x},\lambda)=\big{(}\psi_{m}(\mathbf{x},\lambda)\big{)}_{m=1}^{N_{T}}) for the fluctuating velocity (respectively pressure) associated with each parameter by the method known as snapshots-POD [27]. The POD bases are truncated to $M$ modes for every parameters $\lambda\in\Lambda$ so that the instantaneous fields are approximated as

[TABLE]

In this applcation, $M=10$. We use three methods to interpolate the sets \big{(}\Phi(\mathbf{x},\lambda)\big{)}_{\lambda\in\Lambda} and \big{(}\Psi(\mathbf{x},\lambda)\big{)}_{\lambda\in\Lambda} at the parameter $\lambda^{\star}$: (i) the grassmannian Neville-Aitken interpolation method proposed in subsection 3.2 (labeled Neville), (ii) the method proposed by Amsallem and Farhat in [6] (labeled Amsallem) and (iii)  a naive piecewise affine interpolation of the components associated with the sampled POD bases seen as regular matrices (labeled Standard). We firstly examine the projection error for the fluctuating velocity defined for a given POD basis $\Phi$ as follows:

[TABLE]

where the column of the matrix $\widetilde{\mathbf{V}}_{\lambda}$ are the snapshots of the discrete fluctuating velocity field obtained by the simulation of the HDM for the parameter $\lambda$, $P_{\Phi}\widetilde{\mathbf{V}}=\sum_{m=1}^{M}(\phi_{m}|\widetilde{\mathbf{V}})_{L^{2}(\Omega)}\,\phi_{m}$ denotes the matrix whose columns are the columns of $\widetilde{\mathbf{V}}$ projected onto the subspace engendered by the POD basis $\Phi$, and $\left\|\,\cdot\,\right\|_{F}$ denotes the Frobenius norm. The results are shown in table 2 where the benefit of both interpolations based on the Grassmann manifold is evident when compared with the standard interpolation of the bases components. Also, the proposed grassmannian Neville-Aitken interpolation yields globally a lower error than the method proposed in [6].

Then we construct so-called POD-ROMs by a standard Galerkin projection of the momentum equation onto (i) the POD basis for the velocity and (ii) the gradient of the POD basis for the pressure (see e.g. [28]), with the bases labeled Reference (direct POD of the HDM solution for the target parameter), Neville, Amsallem and Standard. This yields dynamical systems for the temporal coefficients $(a_{m})_{m=1}^{M}$ and $(b_{m})_{m=1}^{M}$ in (15–16). We examine the error associated with the simulation of these dynamical systems defined as follows:

[TABLE]

where the matrix $\widehat{\mathbf{V}}_{\Phi}$ is the reconstruction of the fluctuating velocity $\widehat{\mathbf{v}}-\overline{\mathbf{v}}$ in (15) for the coefficients $(a_{m})_{m=1}^{M}$ obtained by the simulation of the POD-ROM associated with the basis $\Phi$. The results are shown in table 3 where we see that the POD-ROMs built from bases interpolated by the proposed grassmannian Neville-Aitken interpolator yields better results compared with the method proposed in [6].

Finally, we compare the aerodynamic efforts imposed by the fluid on the body, defined as follows:

[TABLE]

where $\Gamma_{\mathrm{disk}}$ denotes the disk boundary, $\sigma$ is the fluid stress tensor, $\mathbf{n}$ is the outward unit vector normal to $\Gamma_{\mathrm{disk}}$, and $\mathrm{d}\ell$ is a line element. Results are shown in figure 2 for the case 2, where we see that the efforts obtained by the POD-ROM built from bases interpolated by the proposed grassmannian Neville-Aitken interpolator are closer to the reference.

4.2 Lid-driven cavity with inflow

Here, we consider the planar flow of an incompressible newtonian fluid in a lid-driven cavity with an inflow and an outflow (see description in figure 3). We define two Reynolds numbers associated with this configuration. The first is associated with the lid velocity $\mathrm{R}_{1}={\frac{\rho\,v_{\mathrm{lid}}\,\ell}{\eta}}$ where $\rho=1$kg.m*-3* is the fluid density, $v_{\mathrm{lid}}=1$m.s*-1* is the magnitude of the lid velocity $\mathbf{v}|_{\Gamma_{1}}=(1,0)$, $\ell=1$m is the length of the cavity’s sides and $\eta=10^{-3}$kg.m*-1*.s*-1* is the dynamic viscosity so that $\mathrm{R}_{1}=10^{3}$. The second is the Reynolds number associated with the velocity at inflow: $\mathrm{R}_{2}={\frac{\rho\,v_{\mathrm{inflow}}\,\ell}{\eta}}$ with $v_{\mathrm{inflow}}$ the magnitude of inflow velocity $\mathbf{v}|_{\Gamma_{\mathrm{IN}}}=(v_{x},0)$. The dimensionless parameter is the ratio of these two Reynolds numbers

[TABLE]

and controls directly the inflow velocity $v_{x}$.

The HDM is derived exactly as in the previous application (subsection 4.1). The mesh includes here 13656 nodes.

Remark that for $\lambda=0$ ($\mathbf{v}|_{\Gamma_{\mathrm{IN}}}=(0,0)$), this application corresponds to the standard lid-driven cavity and that the main vortex develops at another location for $0<\lambda\leq 1$. The test cases are given in table 4. The test procedure is identical to that in previous section 4.1, with the difference that the snapshots are taken from the initial condition over a transient period of 10s. Then, we compute sets of POD bases for the fluctuating velocity and pressure, and we derive interpolated bases by the proposed method and that proposed in [6]. Every bases are truncated to $M=20$ modes.

We firstly examine the projection error for the fluctuating velocity defined in (17). The results are given in table 5, where we see that the error due to the projection of the velocity on the basis interpolated by the proposed method is globally lower than the error due to the projection on the basis interpolated by the method proposed in [6].

Then, we build POD-ROMs from the reference basis Reference and the interpolated bases Neville and Amsallem as described in previous section 4.1. The simulation of these reduce-order dynamical systems yields the temporal coefficients \big{(}a_{i}(t)\big{)}_{i=1}^{M} and \big{(}b_{i}(t)\big{)}_{i=1}^{M} in (15–16). The dynamical systems errors defined in (18) is shown in table 6 for the four cases described in table 4. We see that the proposed method yields the lowest errors.

Finally, we show the result of a bad choice for the reference point $\lambda_{r}$ in the method from [6] in figure 4. This also illustrates the benefit of the independence of the proposed method from any reference point.

4.3 Flows induced by a rotating body

Here, we consider a two dimensional circular spatial domain $\Omega=\Omega_{\mathrm{S}}\cup\Omega_{\mathrm{F}}$ filled with a rotating ellipsoidal body $\Omega_{\mathrm{S}}$ immersed in an incompressible newtonian fluid $\Omega_{\mathrm{F}}$ (see description in figure 5). The parameter is the Reynolds number defined as

[TABLE]

with the density $\rho=1$ (kg.m*-3*), the dynamic viscosity $\eta=0.01$ (kg.m*-1*.s*-1*), the velocity at the ellipse tips ${v_{\infty}=R\,\dot{\theta}}$ (m.s*-1*) and $\ell=2R$ (m) the ellipse principal diameter. The governing equations are derived by extending the incompressible Navier-Stokes equations to the solid domain by the fictitious domains method. Then, the momentum equation and the continuity equation are solved together by a monolithic formulation on the mixed finite element space known as the mini space (i.e. linear vector Lagrange element enriched with cubic vector bubble element for velocity and piecewise linear element for pressure, see [7] for details). The finite-elements mesh includes $52669$ nodes and is not conforming with the body’s boundaries (see figure 5(b)).

The test cases (i.e. the set of sampling parameters and the target parameter) are given in table 7.

The test procedure is identical to that described in section 4.1, except for the construction of the POD-ROM which has to comply with the multiphase description of the fluid–structure interaction (the interested reader is referred to [10]). As in previous section 4.2, we are concerned here with the transient period (half rotation of the ellipse). We compute sets of POD bases for the fluctuating velocity, and we derive interpolated bases by the proposed method and that proposed in [6]. Every bases are truncated here to $M=40$ modes.

The projection error for the fluctuating velocity defined in (17) are given in table 8. The error associated with the reduced-order dynamical systems (in table 9). Globally, the proposed method yields the lowest errors. Additionally, the temporal coefficients \big{(}a_{i}(t)\big{)}_{i=1}^{M} for the fluctuating velocity are shown in figure 6.

5 Conclusion

In this work, we extend the well-known Neville-Aitken’s interpolation algorithm to the Grassmann manifold for the construction of reduced order methods. The recursive character of the algorithm allows us to perform interpolation on the Grassmann manifold using the geodesic barycenter of two points.

Our method presents several advantages with respect to the algorithms developed in [6, 18]. Indeed, it does not require a reference point as in [6]. Moreover, the proposed method is direct and does not require the resolution of a fixed point problem for the minimization process as in [18] in the framework of the IDW method on Grassmann manifolds. Furthermore, our algorithm is more pertinent both in terms of the accuracy of results and the time computation. The performances of our proposed method are illustrated through three independent CFD applications, namely: the Von Karman vortex shedding street, the lid-driven cavity with inflow and the flow induced by a rotating solid.