# Model order reduction by convex displacement interpolation

Simona Cucchiara, Angelo Iollo, Tommaso Taddei, Haysam Telib

► To cite this version:

Simona Cucchiara, Angelo Iollo, Tommaso Taddei, Haysam Telib. Model order reduction by convex displacement interpolation. Journal of Computational Physics, 2024, 514, pp.113230. 10.1016/j.jcp.2024.113230. hal-04808796

HAL Id: hal-04808796

https://inria.hal.science/hal-04808796v1

Submitted on 11 Feb 2025

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
<figure><img src="image_2.png" /></figure>

# Model order reduction by convex displacement interpolation

Simona Cucchiara<sup>a</sup>, Angelo Iollo<sup>b,c,*</sup>, Tommaso Taddei<sup>b,c</sup>, Haysam Telib<sup>a</sup>

## ARTICLE INFO

MSC:

65N30

41A45

35L02

Keywords:

Model order reduction

Nonlinear approximations

## ABSTRACT

We present a nonlinear interpolation technique for parametric fields that exploits optimal transportation of coherent structures of the solution to achieve accurate performance. The approach generalizes the nonlinear interpolation procedure introduced in [Iollo, Taddei, J. Comput. Phys., 2022] to multi-dimensional parameter domains and to datasets of several snapshots. Given a library of high-fidelity simulations, we rely on a scalar testing function and on a point set registration method to identify coherent structures of the solution field in the form of sorted point clouds. Given a new parameter value, we exploit a regression method to predict the new point cloud; then, we resort to a boundary-aware registration technique to define bijective mappings that deform the new point cloud into the point clouds of the neighboring elements of the dataset, while preserving the boundary of the domain; finally, we define the estimate as a weighted combination of modes obtained by composing the neighboring snapshots with the previously-built mappings. We present several numerical examples for compressible and incompressible, viscous and inviscid flows to demonstrate the accuracy of the method. Furthermore, we employ the nonlinear interpolation procedure to augment the dataset of simulations for linear-subspace projection-based model reduction: our data augmentation procedure is designed to reduce offline costs — which are dominated by snapshot generation — of model reduction techniques for nonlinear advection-dominated problems.

## 1. Introduction

### 1.1. Nonlinear interpolation of parametric fields

Despite the many recent contributions, model order reduction (MOR) of parametric problems with compactly-supported features — such as shocks or shear layers — remains an outstanding task for state-of-the-art techniques due to the fundamental inadequacy of linear approximations. The aim of this work is to devise a general — i.e., independent of the underlying equations — interpolation technique for steady-state parametric problems, with emphasis on fluid dynamics applications.

During the past decade, several authors have proposed *Lagrangian methods* to deal with this class of problems [1–8]. We denote by $\mu$ the vector of model parameters in the region $\mathcal{P} \subset \mathbb{R}^p$ and we denote by $\Omega \subset \mathbb{R}^d$ the open computational domain; then, we introduce

<footnote><sup>a</sup> *Optimad Engineering, 10134 Torino, Italy*</footnote>

<footnote><sup>b</sup> IMB, UMR 5251, Univ. Bordeaux, 33400 Talence, France</footnote>

<footnote><sup>c</sup> Inria Bordeaux Sud-Ouest, Team MEMPHIS, 33400 Talence, France</footnote>

<footnote>* Corresponding author at: IMB, UMR 5251, Univ. Bordeaux, 33400 Talence, France.</footnote>

<footnote>E-mail addresses: simona.cucchiara@optimad.it (S. Cucchiara), angelo.iollo@inria.fr (A. Iollo), tommaso.taddei@inria.fr (T. Taddei), haysam.telib@optimad.it (H. Telib).</footnote>
the parametric field of interest $u: \Omega \times \mathcal{P} \to \mathbb{R}^D$ and the solution manifold $m = \{u_\mu := u(\cdot; \mu) : \mu \in \mathcal{P}\}$. Lagrangian approximations rely on the ansatz

$$ \hat{u}_{\mu} = \tilde{u}_{\mu} \circ \Phi_{\mu}^{-1}, \qquad (1a) $$

where $\tilde{u}_\mu$ is a linear (or affine) approximation of the form

$$ \tilde{u}_{\mu}(x) = \sum_{i=1}^{n} \hat{\omega}_{\mu}^{i} \zeta_{i}(x), \quad x \in \Omega, \mu \in \mathcal{P}, \qquad (1b) $$

for proper choices of the weights $\hat{\omega}_{\mu}^1, \dots, \hat{\omega}_{\mu}^n$ and the parameter-independent fields $\zeta_1, \dots, \zeta_n : \Omega \to \mathbb{R}^D$, and $\Phi : \Omega \times \mathcal{P} \to \Omega$ is a suitably-chosen bijection that tracks the coherent structures of the solution; here, $D$ denotes the number of state variables, while $d$ is the spatial dimension. Lagrangian approaches are motivated by the observation that in many problems of interest coherent structures that are troublesome for linear approximations vary smoothly with the parameter and they hence can be tracked through a low-rank parameter-dependent mapping $\Phi$.

In [9], we (Iollo, Taddei) proposed a general method dubbed convex displacement interpolation (CDI) that relies on optimal transportation to perform accurate nonlinear interpolations between solution snapshots; the approach was proposed for databases of two snapshots and one-dimensional parameter domains. Similarly to Lagrangian approaches, CDI relies on the assumption that the location of coherent features of the solution field depends smoothly on the parameter; however, unlike in (1), it does not rely on the definition of a reference configuration where the location of the coherent features is (approximately) frozen. In this work, we discuss an extension of CDI to multi-dimensional parameter domains and to datasets of several snapshots.

## 1.2. Convex displacement interpolation

CDI exploits the standard offline/online decomposition to predict the state map $\mu \in \mathcal{P} \to u_\mu \in m$; Algorithm 1 summarizes the general procedure and highlights the key steps of the methodology that will be discussed in the next sections. Given the solution manifold $m$, we introduce the training set $\mathcal{P}_{\text{train}} = \{\mu^k\}_{k=1}^{n_{\text{train}}} \subset \mathcal{P}$ and the dataset of solutions $\mathcal{D}_{\text{train}} = \{u_\mu : \mu \in \mathcal{P}_{\text{train}}\}$. During the offline (or learning) stage, for each element of $\mathcal{D}_{\text{train}}$, we determine a set of points $X_\mu^{\text{raw}} = \{x_{i,\mu}^{\text{raw}}\}_{i=1}^{N_\mu} \subset \bar{\Omega}$ that describe local features of the solution field that we wish to track: note that the sets $\{X_\mu^{\text{raw}} : \mu \in \mathcal{P}_{\text{train}}\}$ do not necessarily have the same number of elements and are not matched with each other. Then, we define the reference points $X_\mu^{\text{ref}} := \{x_i^{\text{ref}}\}_{i=1}^N$ and the new sensors $X_\mu := \{x_{i,\mu} = x_i^{\text{ref}} + v_{i,\mu}\}_{i=1}^N$ such that the set $X_\mu$ approximates — in a sense to be defined — $X_\mu^{\text{raw}}$ for all $\mu \in \mathcal{P}_{\text{train}}$. During the online stage, given the new parameter value $\mu \notin \mathcal{P}_{\text{train}}$, first, we apply a regression method to predict the new sensor locations $\hat{X}_\mu = \{\hat{x}_{i,\mu} = x_i^{\text{ref}} + \hat{v}_{i,\mu}\}_{i=1}^N$; second, we identify a set of $\kappa$ nearest neighbors $\mathcal{P}_{\text{nn}}^\mu = \{v^i\}_{i=1}^\kappa \subset \mathcal{P}_{\text{train}}$; third, we define the mappings $\{\Phi_v : v \in \mathcal{P}_{\text{nn}}^\mu\}$ such that each mapping $\Phi_v$ is bijective in $\Omega$ and $\Phi_v(\hat{x}_{i,\mu}) \approx \hat{x}_{i,v}$ for $i=1, \dots, N$ and $v \in \mathcal{P}_{\text{nn}}^\mu$; fourth, we return the (generalized) CDI as

$$ \hat{u}_{\mu} = \sum_{v \in \mathcal{P}_{\text{nn}}^{\mu}} \omega_{\mu}^{v} \tilde{u}_{v}, \quad \text{where } \tilde{u}_{v} = u_{v} \circ \Phi_{v}, \qquad (2) $$

for a proper choice of the weights $\{\omega_\mu^v : v \in \mathcal{P}_{\text{nn}}^\mu\}$. Fig. 1 illustrates the computational procedure for a simple one-dimensional problem.

<table><caption>**Algorithm 1:** offline/online decomposition.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td><strong>Offline stage</strong></td><td>performed once</td></tr><tr><td>1: Generate the dataset $\mathcal{D}_{\text{train}} = \{u_\mu : \mu \in \mathcal{P}_{\text{train}}\}$.</td><td></td></tr><tr><td>2: Identify the point clouds $\{X_\mu^{\text{raw}} : \mu \in \mathcal{P}_{\text{train}}\}$.</td><td>section 2.1</td></tr><tr><td>3: Define the template set $X_\mu^{\text{ref}}$ and the sorted point clouds $\{X_\mu : \mu \in \mathcal{P}_{\text{train}}\}$.</td><td>section 2.2</td></tr><tr><td><strong>Online stage</strong></td><td>performed for any $\mu \in \mathcal{P}$</td></tr><tr><td>1: Estimate the new points $\hat{X}_\mu = \{\hat{x}_{i,\mu}\}_{i=1}^N$.</td><td>section 2.3</td></tr><tr><td>2: Select the neighboring parameters $\mathcal{P}_{\text{nn}}^\mu = \{v^i\}_{i=1}^\kappa \subset \mathcal{P}_{\text{train}}$.</td><td>section 2.4</td></tr><tr><td>3: Compute the mappings $\Phi_v$ based on $\hat{X}_\mu$ and $\hat{X}_v$ for all $v \in \mathcal{P}_{\text{nn}}^\mu$.</td><td>section 2.5</td></tr><tr><td>4: Compute the weights $\{\omega_\mu^v : v \in \mathcal{P}_{\text{nn}}^\mu\}$ and return the estimate (2).</td><td>section 2.6</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>Fig. 1. Schematic of the CDI procedure. (a) dataset of snapshots {$u_{\mu^1}, u_{\mu^2}$} and corresponding sensors (shock location), $\hat{x}_{\mu}$ denotes the predicted shock location for the parameter $\mu$. (b) mapped snapshots {$\tilde{u}_{\mu^1}, \tilde{u}_{\mu^2}$} in (2). (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)</figcaption></figure>

The practical implementation of CDI necessitates the careful selection of suitable methods for each phase of both the offline and online stages outlined in Algorithm 1. Below, we present several paradigmatic approaches that serve as a framework, yet they remain open to further refinement and enhancement.

**Sensor selection** There is a broad literature in fluid mechanics and scientific computing on the definition of quantitative criteria for the identification of coherent flow features such as shear layers, shocks, or vortices. Representative examples include shock sensors for high-order methods in CFD [10,11], and local criteria for vortex identification based on swirling strength, and criteria based on the second invariant of the velocity gradient tensor for detecting shear layers [12]. These identification criteria are routinely used in various applications, from aerodynamics to oceanography, allowing the accurate characterization and analysis of complex flow phenomena.

**Point cloud matching** The problem of determining sorted point clouds based on a dataset of unmatched point clouds and a template set is referred to as point-set registration (PSR, [13–17]): PSR methods implicitly introduce a global smooth mapping that coherently deforms the points of the template set $X^{\text{ref}}$ to match — in a sense to be defined — the target cloud $X_{\mu}$, for all $\mu \in \mathcal{P}_{\text{train}}$.

**Regression** The problem of estimating the sensor locations for new parameter values is a standard multi-target regression problem that can be handled through off-the-shelves machine learning routines including radial basis function (RBF, [18]) regression. RBF interpolation/approximation is broadly used in science and engineering due to its ability to cope with sparse datasets and to produce smooth and accurate results. We further notice that RBF approximation is invariant under rototranslations: this feature of RBF approximation is key to prove the frame invariance property of CDI (cf. section 2.7). As in non-intrusive MOR methods (e.g., [19]), we first apply proper orthogonal decomposition (POD, [20,21]) to determine a low-rank representation of the sorted point clouds and then we apply RBF to the POD coefficients.

**Registration** The mapping generated by PSR methods does not fulfill the bijectivity condition $\Phi(\Omega) = \Omega$; in order to avoid extrapolation of full-order solutions outside the domain of definition, it is hence necessary to devise an efficient boundary-aware registration technique that ensures bijectivity in $\Omega$. Note that the registration step should be significantly less expensive than full-order solves in order for the CDI method to be beneficial.

**Choice of the weights** The weights {$\omega_{\mu}^{\nu} : \nu \in \mathcal{P}_{nn}^{\mu}$} should be chosen to ensure relevant properties of the interpolation method (cf. section 2.7). Since our emphasis is on datasets of very modest size, we rely on the inverse distance weighting (IDW, [22]) method.

## 1.3. Contributions and relation to previous works

The past decade has witnessed a surge in the development of nonlinear MOR methods to mitigate the Kolmogorov barrier [23–27]. As discussed in the previous section, our method shares important features with Lagrangian methods based on the ansatz (1); the ansatz (2) is also formally equivalent to front tracking methods proposed in [28,29]. Note also that our approach relies on a piecewise-nonlinear approximation, which is motivated by the need to reduce the mesh interpolation and registration costs that scale linearly with $\kappa$ and also by the particular strategy employed to define the weights {$\omega_{\mu}^{\nu} : \nu \in \mathcal{P}_{nn}^{\mu}$}. In this respect, the method is related to previously-proposed localized MOR techniques [30,31].

The CDI procedure outlined in Algorithm 1 generalizes the method in [9]: in more detail, it reduces to the two-field CDI estimate for one-dimensional parameters and $\kappa = 2$, and a proper choice of the mappings {$\Phi_{\nu} : \nu \in \mathcal{P}_{nn}^{\mu}$} (cf. Remark 1). We observe that, similarly to Lagrangian approximations, CDI relies on geometric mappings to track local features of the solution field and ultimately improve state predictions; on the other hand, CDI does not rely on the definition of a reference domain and the estimate is inherently Eulerian. We note that the definition of the template $X^{\text{ref}}$ is instrumental to the construction of the sorted point clouds {$X_{\mu} : \mu \in \mathcal{P}_{\text{train}}$}, which enable the application of standard regression algorithms during the online stage for the prediction of $\hat{X}_{\mu}$. We also observe that
during the online stage each mapping $\Phi_v$ is computed based on the point clouds $\hat{X}_\mu$ and $\hat{X}_v$ (in lieu of $X_v$): thanks to this choice, we prove that CDI is an interpolation method (cf. section 2.7).

In this work, we rely on well-established techniques for sensor selection, point cloud matching, multi-target regression, and weight selection; on the other hand, we emphasize that the development of boundary-aware registration techniques to determine bijective mappings in the domain $\Omega$ of interest is an active research area with many recent contributions. In this work, we propose a new elasticity-based method that borrows ideas from penalization methods for PDE discretization on Cartesian meshes (see, e.g., [32,33]); we also note that the method shares important features with morphing techniques for mesh adaptation [34,35] and image processing [36,37]. We further consider the optimization-based method first proposed in [8] and further developed in several subsequent works.

The ultimate goal of this work is to devise a nonlinear interpolation strategy that can cope with datasets of very modest size. The evaluation of the CDI estimate (2) is not necessarily online-efficient — that is, the prediction cost scales with the cost of the underlying high-fidelity (HF) mesh — but it is significantly less expensive than the solution to the HF model. In this respect, CDI might be employed to initialize nonlinear solvers for nonlinear problems [38], or to enrich the snapshot dataset that can later be used to train a projection-based reduced-order model (PROM) — in effect, data augmentation [39]. We discuss the application of CDI to data augmentation through the vehicle of a parametric two-dimensional inviscid flow problem.

The outline of the paper is as follows. In section 2, we present the methods for each step of the offline and online stages of Algorithm 1, and we show relevant properties of the interpolation procedure. In section 3, we present a motivating example that illustrates the many features of the approach in a simplified setting. In section 4 we present extensive numerical results for three different test cases that cover compressible and incompressible, viscous and inviscid flows. In section 5, we illustrate the application to data augmentation. Section 6 completes the paper.

## 2. Methodology

### 2.1. Sensor (or feature) selection

We rely on the definition of a problem-dependent scalar testing function $\mathcal{T}$ to identify unmatched point clouds that can be exploited to derive out-of-sample predictions (cf. section 2.3). Given the domain $\Omega$, we define a discrete set of points $P_{hf} \subset \bar{\Omega}$ and we compute $\{\mathcal{T}(x, u_\mu^{hf}) : x \in P_{hf}\}$ for all $\mu \in \mathcal{P}_{train}$, where $u_\mu^{hf}$ denotes the high-fidelity (HF) estimate of the solution field for $\mu \in \mathcal{P}$; then, we set

$$X_{\mu}^{\text{raw}} = \{x_{i,\mu}^{\text{raw}}\} = \{x \in P_{hf} : \mathcal{T}(x, u_{\mu}^{\text{hf}}) \ge t_{\mu}\} \quad (3)$$

for some threshold $t_\mu$ that is specified below.

The choice of the scalar testing function is highly problem-dependent and requires a deep understanding of the underlying physical model and of the features we wish to track. For compressible flow problems in the transonic or supersonic regime, we rely on the Ducros sensor [40,10] to identify the shocks: we denote by $P_{hf} = \{x_k^{hf}\}_{k=1}^{N_e}$ the cell centers of the HF mesh and we compute

$$\mathcal{T}(x_k^{hf}; u_\mu^{hf}) = \sup_{x \in D_k} |\phi(x; u_\mu^{hf})| \quad (4a)$$

where $D_k$ denotes the k-th cell of the mesh with center $x_k^{hf}$, $\phi$ is the Ducros sensor

$$\phi(x; u) := \frac{(-\nabla \cdot \nu)^+}{\sqrt{(\nabla \cdot \nu)^2 + \|\nabla \times \nu\|_2^2 + a^2}} \frac{\|\nabla p\|_2}{p + \epsilon} \|v\|_2, \quad (4b)$$

$v$ is the velocity field, $p$ is the pressure, $a$ is the sound speed and $\epsilon > 0$ is a user-defined tolerance that is set equal to 0.01. Finally, the threshold $t_\mu$ in (3) is chosen equal to the $\gamma_{thr}$ quantile over the training set:

$$t_{\mu} = \text{quantile} \left( \{ \mathcal{T}(x, u_{\mu}^{\text{hf}}) : x \in P_{\text{hf}} \}, \gamma_{\text{thr}} \right). \quad (4c)$$

For the incompressible flow problem of section 4, we rely on the isolines of the streamfunction to identify the recirculation area, for two-dimensional channel flows. In more detail, we define the streamline function

$$\Psi_{\mu}(x = [x_1, x_2]) = \int_{y_{\text{btm}}(x_1)}^{x_2} (u_{\mu}(x_1, s))_1 ds, \quad (5)$$

where $y_{\text{btm}}(x_1) = \inf\{x_2 : [x_1, x_2] \in \Omega\}$ denotes the bottom boundary of the channel; then, we set $\mathcal{T}(x, u_\mu^{hf}) = -\Psi_\mu(x)$ and $t_\mu \equiv 0$; that is,

$$X_{\mu}^{\text{raw}} = \{x \in P_{\text{hf}} : \Psi_{\mu}(x) \le 0\}.$$

We remark that we also tried to exploit the Q-criterion and the triple decomposition of relative motion (TDM) [41]: however, in our experience, these two indicators were not able to distinguish between pure shearing motions and the actual swirling motion of the vortex. The choice of sensor is ultimately based on the modeling approach and the specific requirements of the application.
## 2.2. Point cloud matching

The point clouds in (3) are not necessarily of the same size and are not sorted; therefore, we resort to PSR to find matched point clouds:

$$X_{\mu} = \text{PSR}(X^{\text{ref}}, X_{\mu}^{\text{raw}}), \quad X^{\text{ref}} := X_{\mu}^{\text{raw}}, \quad \mu^{\text{ref}} \in \mathcal{P}_{\text{train}}. \qquad (6)$$

Given the point clouds $X = \{x_i\}_{i=1}^N$ and $Y = \{y_j\}_{j=1}^M$, the problem of PSR consists in finding a map $T: \mathbb{R}^d \to \mathbb{R}^d$ that (approximately) minimizes the distance

$$\text{dist}(Y, T(X)) = \max_{y \in Y} \left( \min_{x \in X} \|y - T(x)\|_2 \right)$$

in a suitable model class; the output of the algorithm is the deformed set of points $T(X)$. A broad range of PSR methods — including the approach considered in this work — relies on a probabilistic interpretation of the problem: these methods rely on the assumption that $\{x_i\}_{i=1}^N$ and $\{y_j\}_{j=1}^M$ are independent identically distributed (iid) samples from a given probability distribution with unknown parameters. As discussed in the introduction, Algorithm 1 can cope with any PSR technique that can be written in the form (6); we refer to a future work for the development of a specialized PSR procedure for MOR. In the remainder of this section, we review the PSR algorithm employed in the numerical experiments of sections 4 and 5.

Following [9] We rely on Gaussian-based PSR. Given the point cloud $X = \{x_i\}_{i=1}^N$, we resort to maximum likelihood estimation (MLE, see, e.g., [42, Chapter 8]) to estimate mean and covariance matrix:

$$\mu_X = \frac{1}{N} \sum_{i=1}^{N} x_i, \quad \Sigma_X = \frac{1}{N} \sum_{i=1}^{N} (x_i - \mu_X)(x_i - \mu_X)^{\top}, \qquad (7)$$

where $(\cdot)^{\top}$ is the transpose operator, and we define the Gaussian model $g_X = \eta(\mu_X, \Sigma_X)$. Similarly, given the point cloud $Y = \{y_j\}_{j=1}^M$, we define the MLE estimates $\mu_Y, \Sigma_Y$ and the corresponding Gaussian model $g_Y = \eta(\mu_Y, \Sigma_Y)$. Finally, we apply the optimal transport map (cf. [43]) from $g_X$ to $g_Y$

$$T_{X,Y}(x) = \mu_Y + \Sigma_X^{-1/2} \left( \Sigma_X^{1/2} \Sigma_Y \Sigma_X^{1/2} \right)^{1/2} \Sigma_X^{-1/2} (x - \mu_X), \qquad (8)$$

to define the deformed point cloud $\tilde{X} = \{\tilde{x}_i = T_{X,Y}(x_i)\}_{i=1}^N$.

## 2.3. Regression of sorted point clouds

We rely on RBF regression to determine the approximation $\mu \in \mathcal{P} \mapsto \hat{X}_{\mu} \in \mathbb{R}^{N \times d}$. Given the sorted point clouds $\{X_{\mu} : \mu \in \mathcal{P}_{\text{train}}\}$, we first apply POD to determine the equivalent representation

$$X_{\mu}^{k} = \sum_{i=1}^{M} Z_{i} \beta_{k}^{i}, \quad \text{with } Z_{1}, \dots, Z_{M} \in \mathbb{R}^{N \times d}, \beta_{1}^{1}, \dots, \beta_{n_{\text{train}}}^{M} \in \mathbb{R},$$

and $M = \min\{N \cdot d, n_{\text{train}}\}$. Then, we apply RBF regression to each coefficient of the expansion separately,

$$\hat{\beta}^{i} = \arg \min_{\beta \in \mathcal{B}_{\phi}} \lambda \|\beta\|_{\mathcal{B}_{\phi}}^{2} + \sum_{k=1}^{n_{\text{train}}} (\beta(\mu^{k}) - \beta_{k}^{i})^{2}, \quad i = 1, \dots, M. \qquad (9)$$

Here, $\mathcal{B}_{\phi}$ denotes the native space associated with the kernel $\phi$ and $\lambda > 0$ is a regularization coefficient. We rely on the Gaussian kernel, and we exploit cross-validation to select $\lambda$ and the kernel width; furthermore, we estimate the out-of-sample R-squared indicator and we keep the coefficients of the expansion for which the latter exceeds the threshold 0.5. In conclusion, we obtain the estimate

$$\hat{X}_{\mu} = \sum_{i=1}^{m} Z_{i} \hat{\beta}_{\mu}^{i}, \qquad (10)$$

with $m \le M$.

## 2.4. Choice of the nearest neighbors

Given the parameter value $\mu \in \mathcal{P}$, we select the parameters $v^1, \dots, v^k \in \mathcal{P}_{\text{train}}$ that minimize the Euclidean distance $\text{dist}(\mu, v) = \|\mu - v\|_2$. The Euclidean distance is invariant under rototranslations of the parameters: as discussed in section 2.7, this property is key to prove the frame indifference of the CDI procedure.
## 2.6. Choice of the weights

We determine the weights {$\omega_{\mu}^{\nu}$ : $v \in \mathcal{P}_{nn}^{\mu}$} using the inverse distance weighting (IDW, [22]) method. Given the nearest-neighbor training parameters $v^1, \dots, v^k \in \mathcal{P}_{train}$ and the new parameter $\mu \in \mathcal{P}$, we consider

$$ \left\{ \begin{array}{ll} \omega_{\mu}^{\nu} = \frac{\bar{\omega}_{\mu}^{\nu}}{\sum_{v' \in \mathcal{P}_{nn}^{\mu}} \bar{\omega}_{v'}^{\nu}}, & \text{with } \bar{\omega}_{\mu}^{\nu} = \frac{1}{\text{dist}^{p}(\mu, v')}, \quad \text{if } \mu \notin \mathcal{P}_{train}; \\ \\ \omega_{\mu}^{\nu} = \begin{cases} 1 & \text{if } \mu = v \\ 0 & \text{otherwise} \end{cases} & \text{if } \mu \in \mathcal{P}_{train}. \end{array} \right. \qquad (14a) $$

The algorithm can cope with arbitrary choices of the distance dist : $\mathcal{P} \times \mathcal{P} \to \mathbb{R}_+$. We here consider the Euclidean distance

$$ \text{dist}(\mu, v) = \|\mu - v\|_2. \qquad (14b) $$

The power parameter $p \ge 1$ influences the relative importance of neighboring parameters: large values of $p$ assign greater influence to values closest to the interpolated parameter. Since we consider relatively small values of $\kappa$, we expect that the sensitivity with respect to $p$ is limited. In the numerical experiments, we set $p = 1$.

**Remark 1.** Let $\mathcal{P} \subset \mathbb{R}$ and let $\kappa = 2$. Given $\mu \in [v^0, v^1]$ with $v_0, v_1 \in \mathcal{P}_{train}$, we define $s = \frac{\mu - v^0}{v^1 - v^0}$; then, exploiting (14a) with $p = 1$, we obtain

$$ \hat{u}_{\mu} = (1 - s)\tilde{u}_{v^0} + s\tilde{u}_{v^1}, $$

which corresponds to the two-field CDI proposed in [9], for a proper choice of the mappings.

## 2.7. Properties of the interpolation procedure

The CDI (2) with the choice of the weights (14) has four noteworthy properties that are of interest for interpolation of fluid mechanics fields.

1. **Interpolation.** Let $\mu \in \mathcal{P}_{train}$. Then, $\hat{u}_{\mu} = u_{\mu}$.

**Proof.** The proof exploits the discussion of the previous sections. Recalling the choice of $\mathcal{P}_{nn}^{\mu}$ in section 2.4, we find $\mu \in \mathcal{P}_{nn}^{\mu}$. Then, exploiting (14a), we find $\hat{u}_{\mu} = \tilde{u}_{\mu} = u_{\mu} \circ \Phi_{\mu}$. Recalling the correction in section 2.3, we have $\hat{X}_{\mu} = X_{\mu}$ and hence (cf. section 2.5) $\Phi_{\mu} = \text{i.d.}$ □

2. **Maximum principle.** Let $b \in \mathbb{R}^D$; then, $\sup_{x \in \Omega} |b \cdot \hat{u}_{\mu}(x)| \le \sup_{x \in \Omega, v \in \mathcal{P}} |b \cdot u_{v}(x)|$ for all $\mu \in \mathcal{P}$.

**Proof.** We have

$$ \begin{aligned} |b \cdot \hat{u}_{\mu}(x)| & \le \sum_{v \in \mathcal{P}_{nn}^{\mu}} \omega_{\mu}^{\nu} |b \cdot \tilde{u}_{v}(x)| & \le \max_{v \in \mathcal{P}_{nn}^{\mu}} \left( \sup_{x \in \Omega} |b \cdot \tilde{u}_{v}(x)| \right) \\ & \stackrel{(ii)}{=} \max_{v \in \mathcal{P}_{nn}^{\mu}} \left( \sup_{x \in \Omega} |b \cdot u_{v}(x)| \right) \le \sup_{v \in \mathcal{P}, x \in \Omega} |b \cdot u_{v}(x)|. \end{aligned} $$

In (i), we applied the definition of CDI, the triangle inequality and the positivity of the weights {$\omega_{\mu}^{\nu}$}$_v$; in (ii), we used the fact that $\sum_{v \in \mathcal{P}_{nn}^{\mu}} \omega_{\mu}^{\nu} = 1$; in (iii), we used that $\Phi_{v}(\Omega) = \Omega$ for all $v \in \mathcal{P}$. □

3. **Minimum principle.** Let $b \in \mathbb{R}^D$; then, $\inf_{x \in \Omega} b \cdot \hat{u}_{\mu}(x) \ge \inf_{x \in \Omega, v \in \mathcal{P}} b \cdot u_{v}(x)$ for all $\mu \in \mathcal{P}$.

**Proof.** Given $x \in \Omega$, $\mu \in \mathcal{P}$, we have indeed that $\omega_{\mu}^{\nu} \ge 0$, $\sum_{v \in \mathcal{P}_{nn}^{\mu}} \omega_{\mu}^{\nu} = 1$, and $b \cdot \tilde{u}_{v}(x) \ge c := \inf_{x \in \Omega, v \in \mathcal{P}} (b \cdot u_{v}(x))$ for all $v \in \mathcal{P}_{nn}^{\mu}$. Therefore,

$$ b \cdot \hat{u}_{\mu}(x) = \sum_{v \in \mathcal{P}_{nn}^{\mu}} \omega_{\mu}^{\nu} (b \cdot \tilde{u}_{v}(x)) \ge c \sum_{v \in \mathcal{P}_{nn}^{\mu}} \omega_{\mu}^{\nu} = c. \quad \square $$

4. **Frame indifference.** Given the parametric field $\mu \mapsto u_{\mu}$, consider the rototranslation $p(\mu) = R\mu + b$ for some rotation matrix $R$ and vector $b \in \mathbb{R}^P$, and the parametric field $p \mapsto v_p = u_{\mu(p)}$. Then, $\hat{u}_{p(\mu)} = \hat{u}_{\mu}$ for all $\mu \in \mathcal{P}$.
**Proof.** It suffices to check that all the steps of the CDI procedure are invariant under roto-translations. First, we observe that RBF interpolation/regression is invariant under rototranslations (cf. [18]). Second, since the choice of the neighboring elements and the choice of the CDI weights (14) are both based on the Euclidean distance, we find that both steps also invariant under roto-translations. Finally, the registration procedure is independent of the parameterization considered and in particular is invariant under roto-translations. □

The maximum and minimum principles can be used to ensure that the CDI prediction is physically meaningful for out-of-sample parameters. To provide a concrete reference, if we apply CDI to the solution to the Euler or Navier-Stokes equations in primitive variables, we can guarantee that density and pressure are strictly positive and do not exceed the maximum value attained in the training set. Similarly, since the mappings {$\Phi_v$}$_v$ in (2) are bijective in $\Omega$, if the far-field conditions are parameter-independent and constant on $\partial\Omega$, we can exploit the same argument to prove that the CDI is consistent with the far-field conditions.

# 3. Motivating example

## 3.1. Model problem

We consider the parametric problem:

$$ egin{cases} -\frac{\partial^2 u_\mu}{\partial x^2} = f_\mu & \text{in } \Omega = (-1, 1), \\ u_\mu(-1) = u_\mu(1) = 0; \end{cases} \quad f_\mu(x) = \frac{1}{\sigma} \exp\left(\frac{(x-\mu)^2}{\sigma^2}\right), \quad (15) $$

with $\mu \in \mathcal{P} = [-0.9, 0.9]$ and $\sigma > 0$. We introduce the space $u = H_0^1(\Omega)$ endowed with the inner product

$$ (w, v) = \int_{-1}^{1} \partial_x w \partial_x v \, dx + wv \, dx, \quad (16a) $$

and the induced norm $|| \cdot || = \sqrt{\langle \cdot, \cdot \rangle}$. Exploiting integration by part, we find the variational problem associated with (15):

$$ \text{find } u_\mu \in \mathcal{l} : a(u_\mu, v) = F_\mu(v) \quad \forall v \in \mathcal{l}, \quad (16b) $$

where

$$ a(w, v) = \int_{-1}^{1} \partial_x w \partial_x v \, dx, \quad F_\mu(v) = \int_{-1}^{1} f_\mu v \, dx. \quad (16c) $$

By tedious but straightforward calculations, we can show that

$$ \frac{\pi}{2+\pi} ||u||^2 \le a(u, u) \le ||u||^2 \quad \forall u \in \mathcal{l}. \quad (17) $$

## 3.2. Linear method

We denote by $\mathcal{Z} \subset \mathcal{l}$ an $n$-dimensional linear subspace of $\mathcal{l}$; we further introduce the Galerkin reduced-order model:

$$ \text{find } \hat{u}_\mu \in \mathcal{Z} : a(\hat{u}_\mu, v) = F_\mu(v) \quad \forall v \in \mathcal{Z}. \quad (18) $$

Exploiting Cèa Lemma for symmetric coercive problems (e.g., [46, Lemma 4.3]) and the estimate (17), we find

$$ \| \hat{u}_{\mu} - u_{\mu} \| \le \sqrt{1 + \frac{2}{\pi} \min_{v \in \mathcal{Z}} \|v - u_{\mu}\|}, \quad \forall \mu \in \mathcal{P}. \quad (19) $$

We conclude that for this model problem Galerkin projection achieves near optimal performance over the entire parameter domain, for any choice of the reduced space $\mathcal{Z}$.

We generate $n_{train} = 15$ snapshots for equispaced parameter values $\mathcal{P}_{train} \subset \mathcal{P}$ and for two values of $\sigma$, $\sigma = 10^{-1}$ and $\sigma = 10^{-3}$. For each value of $\sigma$, we assess performance of the POD Galerkin ROM (18) for three values of the reduced basis size $n$, $n = 5, 10, 15$. Fig. 2(a) shows the behavior of the relative $H^1$ error over $\mathcal{P}$ for three choices of $n$ for $\sigma = 10^{-1}$, while Fig. 2(b) shows the solution and the Galerkin approximations for the parameter that maximizes the prediction error for $n = 15$. Fig. 3 replicates the results for $\sigma = 10^{-3}$.

We observe that the ROM for $n = n_{train}$ provides exact reconstructions for in-sample parameters: this immediately follows from (19). We note, however, that for $\sigma = 10^{-3}$ the error is large — in the order of 20% — and does not decrease as we increase $n$ from $n = 10$ to $n = 15$. We hence conclude that the available snapshot set is not sufficient to achieve an accurate linear representation of the elements of the solution manifold.
<figure><img src="image_4.png" /><figcaption>Fig. 2. Motivating example; linear ROM for $\sigma = 10^{-1}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the error for $n = n_{train}$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 3. Motivating example; linear ROM for $\sigma = 10^{-3}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the error for $n = n_{train}$.</figcaption></figure>

We draw two important conclusions. First, the inadequacy of linear approximations primarily depends on the geometry of the solution manifold — which depends on the PDE model and on the parameterization — rather than on the class of PDEs of interest (hyperbolic, elliptic, parabolic); on the other hand, the PDE model affects the sub-optimality (cf. (19)) of Galerkin projection. Second, we observe that several nonlinear MOR methods (e.g., [24,25]) seek state estimates that are nonlinear functions of the generalized coordinates but, nevertheless, belong to the span of the training snapshots: these methods hence do not overcome the issue of the inaccuracy of the snapshot set.

### 3.3. Convex displacement interpolation

Figs. 4 and 5 show the performance of CDI for the model problem (16) for $\sigma = 10^{-1}$ and $\sigma = 10^{-3}$: Figs. 4(a) and 5(a) show the relative $H^1$ error, while Figs. 4(b) and 5(b) show the behavior of the solution and its estimate for the parameter that maximizes the prediction error. We here set $X_\mu = \{\mu\}$, we consider $\kappa = 2$ and we define the bijective maps $\Phi_\nu$ for $\nu \in \mathcal{P}_{nn}^\mu$ so that $\Phi_\nu(\Omega) = \Omega$, $\Phi_\nu$ is piecewise linear in the intervals $(-1, \mu)$ and $(\mu, 1)$ and $\Phi_\nu(\mu) = \nu$. We notice that the error is below 0.5% for both values of $\sigma$.

### 3.4. Data augmentation

As final test, we perform Galerkin projection based on the parameter-dependent reduced space

$$ \mathcal{Z}_\mu = \operatorname{span} \{ \hat{u}_\mu^{\text{cdi}}, u_{\nu1}, u_{\nu2} \}, \quad \text{with } \mathcal{P}_{nn}^\mu = \{ \nu^1, \nu^2 \}. $$

Figs. 6 and 7 show the results: we observe that the error is roughly three times lower than the error associated with CDI. This result suggests the potential of data augmentation for Galerkin ROMs: nonlinear interpolation is employed to rapidly enhance the accuracy of the reduced space, while Galerkin projection is employed to ensure stable and near-optimal estimates. We further investigate the combination of CDI and data augmentation in the next section.
<figure><img src="image_6.png" /><figcaption>Fig. 4. Motivating example; CDI for $\sigma = 10^{-1}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the prediction error.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 5. Motivating example; CDI for $\sigma = 10^{-3}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the prediction error.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Motivating example; Galerkin ROM with data augmentation for $\sigma = 10^{-1}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the prediction error.</figcaption></figure>

## 4. Nonlinear interpolation of parametric flow fields

In this section, we illustrate the performance of our nonlinear interpolation procedure for three model problems. We rely on the elasticity-based approach for registration (cf. section 2.5.1): in all our tests, we set $\epsilon = 10^{-8}$, $\delta = 50$, and $\eta = 10^{-2}$ in (11a) and (11c); we further consider the time step $\Delta t = 5 \cdot 10^{-3}$ for the integration of (12). We study the sensitivity of the interpolation parameter $s$, for datasets with two parameters $\mu^0, \mu^1$ and $P = 1$. Given the point clouds $X_i = \{x_j^i\}_{j=1}^N \subset \Omega$ for $i = 0, 1$, we define $X(s) = \{(1-s)x_j^0 + sx_j^1\}_{j=1}^N$ and the mappings $\Phi_i(s)$ such that $\Phi_i(X(s), s) \approx X_i$ for $i = 0, 1$; finally, we define the two-field CDI

$$ \hat{u}(s) = (1-s)u_0 \circ \Phi_0(s) + su_1 \circ \Phi_1(s), \qquad (20) $$

and the linear counterpart (convex interpolation, CI)
<figure><img src="image_9.png" /><figcaption>Fig. 7. Motivating example; Galerkin ROM with data augmentation for $\sigma = 10^{-3}$. (a) $H^1$ relative error. (b) comparison of HF and ROM solutions for the parameter that maximizes the prediction error.</figcaption></figure>

$$ \hat{u}(s) = (1 - s)u_0 + su_1. \qquad (21) $$

In the absence of prior information, we here set $s = \frac{\mu - \mu_0}{\mu_1 - \mu_0}$ (“linear relation”, cf. Remark 1); however, we also investigate the behavior of the optimal choice of $s$ as a function of $\mu$,

$$ s_{\mu}^{\text{opt}} = \arg \min_{s \in [0,1]} \|\hat{u}(s) - u_{\mu}\|_{\star}, \qquad (22) $$

where the norm $\| \cdot \|_{\star}$ is introduced below for each model problem.

## 4.1. Compressible flow past a RAE2822 airfoil

### 4.1.1. Model problem

We study the compressible viscous flow past a RAE2822 airfoil for varying inflow Mach number $Ma_{\infty}$ and angle of attack $\alpha$. We model the flow using the compressible Reynolds-averaged Navier-Stokes equations with Menter's SST closure model [47]; we rely on the open-source solver SU2 [48] for snapshot generation. In the remainder, we apply the nonlinear interpolation procedure to estimate the pressure coefficient

$$ C_p = \frac{p - p_{\infty}}{\frac{1}{2} \rho_{\infty} \|v_{\infty}\|^2_{2}}, \qquad (23) $$

where $p, \rho, v$ denote flow pressure, density and velocity, and the subscript $(\cdot)_{\infty}$ denotes the free-stream value of the flow quantity $\cdot$.

### 4.1.2. One-dimensional parameterization

We first study variations of the free-stream Mach number in the range $[0.81, 0.85]$ for $\alpha = 2^{\circ}$. Fig. 8(a) shows the computational domain; Figs. 8(b) and (c) illustrate the behavior of the pressure coefficients for $Ma_{\infty} = 0.81$ and $Ma_{\infty} = 0.85$. We observe that the solution develops two shocks on the lower and upper sides of the airfoil whose position and intensity strongly depend on the parameter value.

We rely on the Ducros sensor (4) to identify the markers, with $\gamma_{thr} = 99.6\%$ (cf. (4c)). This sensor is designed to detect regions of strong compression with abrupt changes in flow properties, such as density, pressure, or velocity gradients, which are characteristic of shock waves. We note that the numerical results depend on the value of $\gamma_{thr}$: excessively large values of $\gamma_{thr}$ prevent the sensor from properly capturing the two shock waves, while excessively small values of $\gamma_{thr}$ lead to the selection of cells where shocks are not present. Nevertheless, in our experiments, we found that the optimal choice of $\gamma_{thr}$ is nearly insensitive to the choice of the parameters: we can indeed consider the same $\gamma_{thr}$ for all snapshots in the training set.

Fig. 9 shows the marked cells of the mesh for two parameter values. Since the flow exhibits two distinct coherent structures (shocks), we apply the procedure of section 2.2 to the lower and upper sides of the domain separately (see also [9, section 5.2.2]).

Fig. 10 investigates the performance of CDI (20) and CI (21) over $Ma_{\infty} \in [0.81, 0.85]$ for linear choice of the interpolation parameter $s$ on the boundary $\Gamma_2$. We observe that CDI preserves the structure of the solution and is also significantly more accurate. Fig. 11 compares the behavior of the optimal choice of $s$ (“experimental”) with the linear relationship — optimality is measured with respect to the $L^2(\Gamma_2)$ norm over the airfoil (“integral error”). We observe that for this model problem $s_{\mu}^{\text{opt}}$ varies linearly with respect to the parameter $\mu = Ma_{\infty}$.

### 4.1.3. Two-dimensional parameterization

We consider the problem of estimating the solution for $\mu = (Ma_{\infty}, \alpha) = [0.83, 2]$, based on the training points $\mu^1 = [0.82, 2]$, $\mu^2 = [0.84, 1.5]$ and $\mu^3 = [0.84, 2.5]$. We compare the performance of convex interpolation and CDI. Given the (very) modest amount of datapoints, we rely on linear interpolation for the regression step (cf. Algorithm 1, Line 1, Online stage). Fig. 12 compares the
<figure><img src="image_10.png" /><figcaption>Fig. 8. Compressible flow past a RAE2822 airfoil. (a) computational domain. (b)-(c) pressure coefficients for Ma∞ = 0.81 and Ma∞ = 0.85, α = 2°.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 9. Compressible flow past a RAE2822 airfoil. Marked cells for Ma∞ = 0.81 and Ma∞ = 0.85, α = 2°.</figcaption></figure>

predictions of truth and predicted pressure coefficient for the out-of-sample parameter μ. We observe that nonlinear interpolation is qualitatively more accurate, especially in the neighborhood of the shock. We observe that the L²(Ω) error is 5.62% for CI and 4.10% for CDI; the L²(Γ₂) error is 2.69% for CI and 2.14% for CDI. We emphasize that the L² metric does not fully capture local improvements in the regions where the shocks develop.

## 4.2. Incompressible flow past a backward facing step

### 4.2.1. Model problem

We study the incompressible flow past a backflow facing step for varying inflow Reynolds number Re and step's inclination angle α; we model the flow using the incompressible Navier-Stokes equations and we rely on the open-source software OpenFOAM [49].
<figure><img src="image_12.png" /></figure>

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

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

<figure><img src="image_15.png" /><figcaption>Fig. 10. Compressible flow past a RAE2822 airfoil, one-dimensional parameterization. Behavior of truth and predicted pressure coefficient on lower (LOW) and upper (UP) sides of the blade, for two values of Ma∞ and α = 2°. (a)-(b) convex interpolation. (c)-(d) convex displacement interpolation.</figcaption></figure>

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

<figure><img src="image_17.png" /><figcaption>Fig. 11. Compressible flow past a RAE2822 airfoil, one-dimensional parameterization. (a) behavior of the optimal value of s in (22) for CDI and CI. (b) behavior of the integral error over the airfoil.</figcaption></figure>

Mesh morphing is based on the Optimad proprietary software mimic.<fnref n="1" /> The coherent structures of interest are recirculation bubbles (vortices), whose extensions and positions vary as a function of the parameters. We apply the nonlinear interpolation procedure to the velocity field.

<footnote n="1">¹ See https://optimad.github.io/mimmo/.</footnote>
<figure><img src="image_18.png" /><figcaption>Fig. 12. Compressible flow past a RAE2822 airfoil, two-dimensional parameterization. Behavior of truth and predicted pressure coefficient for an out-of-sample configuration. (a) convex interpolation. (b) convex displacement interpolation.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 13. Incompressible flow past a backward facing step. Velocity magnitude for two values of Re and $\alpha = 90^{\circ}$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 14. Incompressible flow past a backward facing step. Marked cells for two values of Re and $\alpha = 90^{\circ}$.</figcaption></figure>

### 4.2.2. One-dimensional parameterization

We first study variations of the inflow Reynolds number in the range [50, 200] for $\alpha = 90^{\circ}$. Fig. 13 shows the velocity magnitude for two parameter values, while Fig. 14 shows the recirculation area as predicted by the indicator (5). We observe that the recirculation area is very sensitive to the value of the inflow Reynolds number.

Fig. 15 compares of CDI (20) and CI (21) for Re = 100. We observe that the nonlinear interpolation is qualitatively more accurate to represent the recirculation area. Fig. 16 investigates the behavior of the optimal choice of $s$ (“experimental”) and compares it with the linear relationship — optimality is measured with respect to the $L^2(\Omega)$ norm. We observe that $s_{\mu}^{\text{opt}}$ depends almost linearly on Re;
<figure><img src="image_21.png" /><figcaption>**Fig. 15.** Incompressible flow past a backward facing step, one-dimensional parameterization. Predictions of velocity streamlines for $Re = 100$ based on the high-fidelity model, nonlinear interpolation and linear interpolation.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 16.** Incompressible flow past a backward facing step, one-dimensional parameterization. (a) behavior of the optimal value of $s$ in (22) for CDI and CI. (b) behavior of the relative $L^2(\Omega)$ error for the velocity field.</figcaption></figure>

we further observe that the improvement due to the nonlinear interpolation is less significant than for the previous case; we again emphasize that the global $L^2$ norm does not properly capture local-in-space features of the error.

### 4.2.3. Two-dimensional parameterization

We consider the problem of estimating the solution for $\mu = (\mathrm{Re}, \alpha) = [120, 70^\circ]$ based on the training points $\mu^1 = [100, 50^\circ]$, $\mu^2 = [100, 90^\circ]$ and $\mu^3 = [180, 70^\circ]$. We compare performance of convex interpolation and CDI in terms of streamlines and horizontal velocity prediction. As in the previous test case, prediction of the markers' locations for out-of-sample configurations is performed through linear interpolation. Fig. 17 shows the results.

We observe that the improvement due to nonlinear interpolation in terms of global $L^2(\Omega)$ relative error in velocity prediction is very modest (1.31% vs 1.29%). However, we notice that the streamlines of the flow field that are predicted through nonlinear interpolation are qualitatively more accurate than the ones associated with linear interpolation; similarly, the prediction of the horizontal velocity based on the nonlinear model over the horizontal slice $\{x: x_1 \in (0.1, 0.7), x_2 = -0.01\}$ is significantly more accurate in the recirculation region, especially in the proximity of the minimum of the horizontal velocity.
<figure><img src="image_23.png" /><figcaption>**Fig. 17.** Incompressible flow past a backward facing step two-dimensional parameterization. Comparison of high-fidelity, nonlinear interpolation, and linear interpolation predictions. (a)-(b)-(c) streamlines of the velocity field. (d) horizontal velocity over the horizontal slice {$x: x_1 \in (0.1, 0.7), x_2 = -0.01$}.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 18.** Compressible flow past a sphere. Pressure coefficients for two values of $Ma_\infty$.</figcaption></figure>

## 4.3. Compressible flow past a sphere

### 4.3.1. Model problem

We study the compressible supersonic inviscid flow past a sphere for inflow Mach number $Ma_\infty \in [2, 5]$. We model the flow using the compressible Euler equations for ideal gases and we rely on SU2 for snapshot generation. As in section 4.1, we apply the nonlinear interpolation strategy to estimate the pressure coefficient (23) and we rely on the Ducros sensor with $\gamma_{thr} = 99.6\%$ to identify the coherent structure. Fig. 18 shows the behavior of $C_p$ for $Ma_\infty = 2$ and $Ma_\infty = 5$: we observe that the solution exhibits a bow shock ahead of the sphere whose location depends on the parameter. Fig. 19 shows the marked cells for the same parameter values.
<figure><img src="image_25.png" /><figcaption>Fig. 19. Compressible flow past a sphere. Marked cells for Ma∞ = 2 and Ma∞ = 5.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Fig. 20. Compressible flow past a sphere. (a) behavior of the optimal value of s in (22) for CDI and CI. (b) behavior of the relative L²(Ω) error for the pressure coefficient.</figcaption></figure>

Fig. 20 illustrates the performance of CDI (20) and CI (21) for Ma∞ ∈ [2, 5]: Fig. 20(a) investigates the behavior of the optimal value of s in (22) for CDI and CI; Fig. 20(b) shows the behavior of the relative L²(Ω) error for the pressure coefficient. Interestingly, we observe that the behavior of $s_{\mu}^{\text{opt}}$ is strongly nonlinear with respect to the Mach number. As for the previous cases, nonlinear interpolation is superior to linear interpolation, especially when the parameter s in (20) is chosen in an optimal manner.

## 5. Data augmentation

### 5.1. Model problem: inviscid flow past a bump

We consider the problem of estimating the solution to the two-dimensional Euler equations in the transonic regime in a channel with a Gaussian bump. We denote by $ρ$ the density, by $v$ the velocity field, by $E$ the total energy, and by $p = (γ-1)(E - \frac{1}{2}ρ||v||_2^2)$ the pressure, where $γ > 0$ is the ratio of specific heats; we further denote by $u = [ρ, ρv, E]$ the vector of conserved variables. We introduce the computational domain $Ω_μ = \{x ∈ (-1.5, 1.5) × (0, 0.8) : x_2 > he^{-25x_1^2}\}$ where $h > 0$ is a given parameter. We consider the conservation law:

$$ \nabla \cdot F(u^{\text{true}}) = 0, \text{ where } F(u) = \begin{bmatrix} \rho v^T \\ \rho u v^T \\ (E+p)v^T \end{bmatrix}, \qquad (24) $$

completed with wall boundary conditions on top and bottom boundaries, subsonic inlet condition (total temperature, total pressure and flow direction) at the left boundary and subsonic outlet condition (static pressure) at the right boundary. We express the free-stream field $u_∞$ in terms of the Mach number Ma_∞,
<figure><img src="image_27.png" /><figcaption>Fig. 21. Inviscid flow over a Gaussian bump. Visualization of the Mach field for two parameter values.</figcaption></figure>

$$T_{\infty} = 1, \quad p_{\infty} = \frac{1}{\gamma}, \quad \rho_{\infty} = 1, \quad u_{\infty} = [Ma_{\infty}, 0],$$

where $T = \frac{p}{R\rho}$ is the temperature and $R = \gamma - 1$ is the (non-dimensional) specific gas constant. Finally, we introduce the parameter vector $\mu = [h, Ma_{\infty}]$ and the parameter region $\mathcal{P} = [0.05, 0.065] \times [0.68, 0.78]$. Figs. 21(a) and (b) show the behavior of the Mach field for two values of the parameters in $\mathcal{P}$.

## 5.2. Geometric parameterization and high-fidelity solver

The computational domain $\Omega_{\mu}$ depends on the geometric parameter $h$; therefore, we should introduce a geometric mapping to recast the problem over a parameter-independent configuration. We here resort to a Gordon-Hall transformation: we denote by $\Phi_{\text{geo},\mu}$ the geometric mapping from the rectangle $\tilde{\Omega} = (-1.5, 1.5) \times (0, 0.8)$ to the physical domain $\Omega_{\mu}$ (see, e.g., [45, section 2]). We further define the domain $\Omega := \Omega_{\tilde{\mu}}$, where $\tilde{\mu}$ is the centroid of $\mathcal{P}$ and the additional geometric mapping $\Psi_{\text{geo},\mu} := \Phi_{\text{geo},\mu} \circ \Phi_{\text{geo},\tilde{\mu}}^{-1} : \Omega \to \Omega_{\tilde{\mu}}$.

The reference configuration $\tilde{\Omega}$ is used below to define the sensors and the regression procedure; on the other hand, the HF mesh is defined in the domain $\Omega$. The latter choice is justified by the observation that for small variations of $h$ the mapping $\Psi_{\text{geo}}$ is (significantly) closer to the identity map than $\Phi_{\text{geo}}$; therefore, defining the mesh in $\Omega$ offers significantly enhanced control over the quality of the deformed mesh.

Simulations are performed in Matlab 2022a [50] based on an in-house code, and executed over a commodity Linux workstation (RAM 32 GB, Intel i7 CPU 3.20 GHz x 12). We consider a P2 mesh with $N_e = 12748$ elements; given a new value of the parameter $\mu$, we define the deformed mesh using the geometric mapping $\Psi_{\text{geo}}$, and then we call the DG solver to compute the solution field. We measure performance and we compute reduced spaces in the reference domain $\Omega$; we denote by $(\cdot, \cdot)$ the $[L^2(\Omega)]^4$ inner product and by $\|\cdot\| = \sqrt{(\cdot, \cdot)}$ the induced norm.

## 5.3. Model order reduction

We rely on a projection-based ROM based on a least-square Petrov-Galerkin (LSPG, [51]) formulation and empirical quadrature (EQ, [52,53]) hyper-reduction; we refer to [45] for a detailed presentation of the method.

## 5.4. Convex displacement interpolation

We apply Algorithm 1 to construct the CDI estimate.

**Feature selection.** We rely on the Ducros sensor (4) to identify the elements of the mesh that contain the shock. We set $\gamma_{\text{thr}} = 99\%$ in (4c) to determine the threshold $t_{\mu}$. In order to take into account geometric variations, we define the “mapped” point clouds $\tilde{X}_{\mu}^{\text{raw}} = \{\Phi_{\text{geo},\mu}^{-1}(x_{i,\mu}^{\text{raw}})\}_{i=1}^{N_{\mu}} \subset \tilde{\Omega}$, where $\Phi_{\text{geo},\mu} : \tilde{\Omega} \to \Omega_{\mu}$ is the Gordon-Hall map. Then, we set $\tilde{X}^{\text{ref}} = \tilde{X}_{\mu}^{\text{raw}}$ where $\mu^{\text{ref}}$ is the parameter in $\mathcal{P}_{\text{train}}$ that minimizes the distance from the centroid of $\mathcal{P}$, and we rely on Gaussian-based PSR to determine the sorted point clouds $\{\tilde{X}_{\mu} : \mu \in \mathcal{P}_{\text{train}}\}$.

**Regression.** We apply the regression procedure described in section 2.3 to find the RBF model $\mu \in \mathcal{P} \mapsto \tilde{X}_{\mu}^{\star} = \{\tilde{x}_{i,\mu}^{\star}\}_{i=1}^{N_{\mu}} \subset \tilde{\Omega}$. Then, we define the point cloud in physical domain using the geometric mapping $\Phi_{\text{geo}}: \tilde{X}_{\mu} = \{\tilde{x}_{i,\mu} = \Phi_{\text{geo},\mu}(\tilde{x}_{i,\mu}^{\star})\}_{i=1}^{N_{\mu}}$. The choice of applying RBF and PSR in the reference (parameter-independent) domain $\tilde{\Omega}$ robustifies the learning process and ultimately improves prediction performance.

**Nearest neighbors, registration, and choice of the weights.** We build the CDI based on $\kappa = 4$ neighbors; we consider the Euclidean distance to identify the neighbors and also in the IDW strategy for the weights (cf. section 2.6). Finally, we resort to optimization-based
<figure><img src="image_28.png" /><figcaption>Fig. 22. Nonlinear interpolation. Behavior of the average relative $L^2$ error for LSPG, CDI and CI for snapshot sets of increasing size $n_{train}$.</figcaption></figure>

registration (cf. (13)) to determine the mapping $\Phi$: we consider the penalty term $\hat{f}_{jac}$ ([54, section 4.3]) that penalizes small values of the Jacobian determinant, and we rely on a piecewise-polynomial search space with $M = 198$ degrees of freedom.

## 5.5. Data augmentation procedure

We denote by $\mathcal{P}_{train} = \{\mu^i\}_{i=1}^{n_{train}} \subset \mathcal{P}$ the training set of parameters for which HF simulations are available and by $\mathcal{P}_{train,cdi} = \{\nu^i\}_{i=1}^{n_{train,cdi}} \subset \mathcal{P}$ an additional set of parameters. We first define the Lagrangian space $\mathcal{Z}_0 = \text{span}\{u_{\mu^i}^{hf}\}_{i=1}^{n_{train}}$ and the projection operator $\Pi_{\mathcal{Z}_0} : [L^2(\Omega)]^4 \to \mathcal{Z}_0$; then, we define the POD space

$$ \mathcal{Z}_n = \mathcal{Z}_0 \oplus \mathcal{Z}_{lf}, \quad \mathcal{Z}_{lf} = \text{POD}(\mathcal{D}_{train,lf}, n - n_{train}, (\cdot, \cdot)), \quad (25) $$

where $\oplus$ indicates the direct sum of orthogonal spaces, $\mathcal{D}_{train,lf} = \{\hat{u}_{\mu} - \Pi_{\mathcal{Z}_0} \hat{u}_{\mu} : \mu \in \mathcal{P}_{train,cdi}\}$ and $\mathcal{Z}_{lf}$ is the POD space of size $n - n_{train}$ associated with the snapshots $\mathcal{D}_{train,lf}$ and the $L^2(\Omega)$ inner product $(\cdot, \cdot)$. We observe that the space $\mathcal{Z}_n$ contains the HF training set: potential inaccuracies of the CDI estimates do not pollute the performance of the reduced space for $\mu \in \mathcal{P}_{train}$.

## 5.6. Numerical results: nonlinear interpolation

We first investigate the performance of CDI for this model problem. Towards this end, we study the behavior of the average relative $L^2(\Omega)$ state prediction error of the projection-based linear ROM (LSPG), CDI, and convex interpolation (CI), for tensorized grids of increasing size $n_{train}$; in all our experiments we set $\kappa = 4$ for CDI and CI, while the reduced space of LSPG includes all the training snapshots ($n = n_{train}$). Fig. 22 shows the results.

We observe that CDI outperforms CI for all parameter values and all training set sizes $n_{train}$. Interestingly, CDI is more accurate than LSPG for $n_{train} = 4$: this can be explained by observing that, unlike the LSPG ROM, CDI provides estimates that do not belong to the span of the dataset of HF simulations; for small values of $n_{train}$, CDI can hence compensate the fact that it does not exploit the knowledge of the PDE model for prediction.

## 5.7. Numerical results: data augmentation

Fig. 23 illustrates the performance for $n_{train} = 9$, $n_{train,cdi} = 20$. We consider a tensorized 3 × 3 grid for $\mathcal{P}_{train}$ and we generate $\mathcal{P}_{train,cdi}$ using latin hyper-cube sampling [55]. We consider a POD space in (25) of dimension $n = 29$; as for the previous experiment, we assess performance based on a dataset of $n_{test} = 20$ out-of-sample parameters. Fig. 23(a) shows boxplots of the relative $L^2(\Omega)$ error for LSPG ROM, CDI and LSPG ROM based on $n = 9$ modes (i.e., without data augmentation) (LSPG-0). We observe that CDI is not as accurate as LSPG-0; however, it can be used to augment the dataset of simulations and ultimately contribute to reduce the prediction error. Fig. 23(b) shows the projection error over the test for the POD space built using (25) (“mixed training”), the POD space based on $n_{train} + n_{train,cdi}$ HF simulations (“HF training”) and the POD space based $n_{train}$ HF simulations (“HF training ($n_{train} = 9$)”: we observe that the mixed training strategy based on CDI is not as accurate as HF training, particularly for large values of $n$; however, it enables significant reduction of the projection error without resorting to new simulations.

Fig. 24 replicates the results of Fig. 23 for $n_{train} = 4$, $n_{train,cdi} = 25$. Note that CDI is more accurate than LSPG-0, which is consistent with the results of Fig. 22.

Table 1 provides an overview of the computational costs. The HF solver is roughly 14 times slower than CDI and 625 times slower than the LSPG ROM; the construction of the LSPG ROM, which comprises the construction of the trial and test spaces and hyper-reduction, is roughly as expensive as a single full-order solve. We remark that approximately 99% of the total cost of CDI is associated with mesh interpolation: we hence envision that more sophisticated implementations of this step (see, e.g., [56]) might lead to a significant reduction of the overall CDI cost.
<figure><img src="image_29.png" /><figcaption>Fig. 23. Data augmentation $n_{train} = 9$, $n_{train,cdi} = 20$. Behavior of LSPG ROM with CDI-based data augmentation. (a) boxplots of relative $L^2$ error for LSPG ROM, CDI and projection error ($n = 29$). (b) projection error.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Fig. 24. Data augmentation $n_{train} = 4$, $n_{train,cdi} = 25$. Behavior of LSPG ROM with CDI-based data augmentation. (a) boxplots of relative $L^2$ error for LSPG ROM, CDI and projection error ($n = 29$). (b) projection error.</figcaption></figure>

<caption>**Table 1**</caption>

Data augmentation. Wall-clock costs of HF and MOR procedures ($n_{train} = 9$, $n_{train,cdi} = 20$).

<table><thead><tr><th></th><th>cost [s]</th></tr></thead><tbody><tr><td>HF solver (avg)</td><td>481.21</td></tr><tr><td>ROM generation (data compression+hyperreduction)</td><td>431.24</td></tr><tr><td>Convex displacement interpolation (avg)</td><td>33.49</td></tr><tr><td>LSPG ROM (n = 29, avg)</td><td>0.77</td></tr></tbody></table>

## 6. Conclusions

We presented a general (i.e., independent of the underlying equations) nonlinear interpolation procedure for parametric fields with local compactly-supported coherent structures; the approach is dubbed convex displacement interpolation (CDI) and extends the previous proposal of [9]. The methodology comprises five main building blocks that were all discussed in section 2: sensor selection, point cloud matching, regression of sorted point clouds, boundary-aware registration, and piecewise-multivariate interpolation. CDI assumes that coherent structures of the solution field — such as shocks, shear layers, vortices — vary smoothly with the parameter. The numerical experiments suggest that this assumption holds for a broad and diverse range of parametric problems in continuum mechanics.

We illustrated the application of CDI to data augmentation. Since CDI generates estimates that do not belong to the span of the training snapshots, we can exploit predictions for out-of-sample parameters to augment the dataset of simulations that are later used for the construction of projection-based ROMs. We envision that combination of lower and higher fidelity snapshots for MOR training has the potential to significantly reduce offline costs, and ultimately enable the application of MOR techniques to complex engineering problems for which extensive parameter explorations based on HF solves are not practical.
CRediT authorship contribution statement

**Simona Cucchiara:** Conceptualization, Methodology, Writing – original draft, Writing – review & editing. **Angelo Iollo:** Conceptualization, Methodology, Writing – original draft, Writing – review & editing. **Tommaso Taddei:** Conceptualization, Methodology, Writing – original draft, Writing – review & editing. **Haysam Telib:** Conceptualization, Methodology, Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to
influence the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgements

The authors acknowledge the support by European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA).

Appendix A. Supplementary material

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jcp.2024.113230.

References

[1] D.S. Ching, P.J. Blonigan, F. Rizzi, J.A. Fike, Model reduction of hypersonic aerodynamics with residual minimization techniques, in: AIAA SCITECH 2022 Forum, 2022, p. 1247.

[2] A. Iollo, D. Lombardi, Advection modes by optimal mass transfer, Phys. Rev. E 89 (2) (2014) 022923.

[3] M.A. Mirhoseini, M.J. Zahr, Model reduction of convection-dominated partial differential equations via optimization-based implicit feature tracking, J. Comput. Phys. 473 (2023) 111739.

[4] R. Mojgani, M. Balajewicz, Arbitrary Lagrangian Eulerian framework for efficient projection-based reduction of convection dominated nonlinear flows, in: APS Division of Fluid Dynamics Meeting Abstracts, 2017, pp. M1–008.

[5] R. Mojgani, M. Balajewicz, Low-rank registration based manifolds for convection-dominated PDEs, in: Proceedings of the AAAI Conference on Artificial Intelligence, vol. 35, 2021, pp. 399–407.

[6] M. Ohlberger, S. Rave, Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing, C. R. Math. 351 (23–24) (2013) 901–906.

[7] T. Taddei, S. Perotto, A. Quarteroni, Reduced basis techniques for nonlinear conservation laws, ESAIM: Math. Model. Numer. Anal. 49 (3) (2015) 787–814.

[8] T. Taddei, A registration method for model order reduction: data compression and geometry reduction, SIAM J. Sci. Comput. 42 (2) (2020) A997–A1027.

[9] A. Iollo, T. Taddei, Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models, J. Comput. Phys. 471 (2022) 111671.

[10] F. Nicoud, F. Ducros, Subgrid-scale stress modelling based on the square of the velocity gradient tensor, Flow Turbul. Combust. 62 (3) (1999) 183–200.

[11] P.-O. Persson, J. Peraire, Sub-cell shock capturing for discontinuous Galerkin methods, in: 44th AIAA Aerospace Sciences Meeting and Exhibit, 2006, p. 112.

[12] J. Jeong, F. Hussain, On the identification of a vortex, J. Appl. Fluid Mech. 285 (1995) 69–94.

[13] R. Horaud, F. Forbes, M. Yguel, G. Dewaele, J. Zhang, Rigid and articulated point registration with expectation conditional maximization, IEEE Trans. Pattern Anal. Mach. Intell. 33 (3) (2010) 587–602.

[14] J. Ma, W. Qiu, J. Zhao, Y. Ma, A.L. Yuille, Z. Tu, Robust L₂ E estimation of transformation for non-rigid registration, IEEE Trans. Signal Process. 63 (5) (2015) 1115–1129.

[15] J. Ma, J. Wu, J. Zhao, J. Jiang, H. Zhou, Q.Z. Sheng, Nonrigid point set registration with robust transformation learning under manifold regularization, IEEE Trans. Neural Netw. Learn. Syst. 30 (12) (2018) 3584–3597.

[16] B. Maiseli, Y. Gu, H. Gao, Recent developments and trends in point set registration methods, J. Vis. Commun. Image Represent. 46 (2017) 95–106.

[17] B. Zitova, J. Flusser, Image registration methods: a survey, Image Vis. Comput. 21 (11) (2003) 977–1000.

[18] H. Wendland, Scattered Data Approximation, vol. 17, Cambridge University Press, 2004.

[19] M. Guo, J.S. Hesthaven, Reduced order modeling for nonlinear structural analysis using Gaussian process regression, Comput. Methods Appl. Mech. Eng. 341 (2018) 807–826.

[20] L. Sirovich, Turbulence and the dynamics of coherent structures. I. coherent structures, Q. Appl. Math. 45 (3) (1987) 561–571.

[21] S. Volkwein, Model reduction using proper orthogonal decomposition, Lecture Notes, Institute of Mathematics and Scientific, Computing, University of Graz, http://www.uni-graz.at/imawww/volkwein/POD.pdf, 2011, 1025.

[22] D. Shepard, A two-dimensional interpolation function for irregularly-spaced data, in: Proceedings of the 1968 23rd ACM National Conference, 1968, pp. 517–524.

[23] D. Amsalem, C. Farhat, Interpolation method for adapting reduced-order models and application to aeroelasticity, AIAA J. 46 (7) (2008) 1803–1813.

[24] J. Barnett, C. Farhat, Quadratic approximation manifold for mitigating the Kolmogorov barrier in nonlinear projection-based model order reduction, J. Comput. Phys. 464 (2022) 111348.

[25] J.L. Barnett, C. Farhat, Y. Maday, Neural-network-augmented projection-based model order reduction for mitigating the Kolmogorov barrier to reducibility of CFD models, arXiv preprint, arXiv:2212.08939, 2022.

[26] K. Lee, K.T. Carlberg, Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders, J. Comput. Phys. 404 (2020) 108973.

[27] B. Peherstorfer, Model reduction for transport-dominated problems via online adaptive bases and adaptive sampling, SIAM J. Sci. Comput. 42 (5) (2020) A2803–A2836.

[28] P. Krah, S. Büchholz, M. Häringer, J. Reiss, Front transport reduction for complex moving fronts, J. Sci. Comput. 96 (28) (2023).
[29] J. Reiss, P. Schulze, J. Sesterhenn, V. Mehrmann, The shifted proper orthogonal decomposition: a mode decomposition for multiple transport phenomena, SIAM J. Sci. Comput. 40 (3) (2018) A1322–A1344.

[30] D. Amsallem, M.J. Zahr, C. Farhat, Nonlinear model order reduction based on local reduced-order bases, Int. J. Numer. Methods Eng. 92 (10) (2012) 891–916.

[31] J.L. Eftang, A.T. Patera, E.M. Rønquist, An “hp” certified reduced basis method for parametrized elliptic partial differential equations, SIAM J. Sci. Comput. 32 (6) (2010) 3170–3200.

[32] F. Chantalat, C.-H. Bruneau, C. Galusinski, A. Iollo, Level-set, penalization and Cartesian meshes: a paradigm for inverse problems and optimal design, J. Comput. Phys. 228 (17) (2009) 6291–6315.

[33] P. Angot, C.-H. Bruneau, P. Fabrie, A penalization method to take into account obstacles in incompressible viscous flows, Numer. Math. 81 (4) (1999) 497–520.

[34] T.E. Tezduyar, M. Behr, S. Mittal, J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain/space-time procedure: II. computation of free-surface flows, two-liquid flows, and flows with drifting cylinders, Comput. Methods Appl. Mech. Eng. 94 (3) (1992) 353–371.

[35] P. Tonon, R.A.K. Sanches, K. Takizawa, T.E. Tezduyar, A linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion, Comput. Mech. 67 (2021) 413–434.

[36] V. Camion, L. Younes, Geodesic interpolating splines, in: International Workshop on Energy Minimization Methods in Computer Vision and Pattern Recognition, Springer, 2001, pp. 513–527.

[37] Y. Cao, M.I. Miller, R.L. Winslow, L. Younes, Large deformation diffeomorphic metric mapping of vector fields, IEEE Trans. Med. Imaging 24 (9) (2005) 1216–1230.

[38] T. Kadeethum, D. O'Malley, F. Ballarin, I. Ang, J.N. Fuhg, N. Bouklas, V.L. Silva, P. Salinas, C.E. Heaney, C.C. Pain, Enhancing high-fidelity nonlinear solver with reduced order model, Sci. Rep. 12 (1) (2022) 20229.

[39] F. Bernard, A. Iollo, S. Riffaud, Reduced-order model for the BGK equation based on POD and optimal transport, J. Comput. Phys. 373 (2018) 545–570.

[40] D. Modesti, S. Pirozzoli, A low-dissipative solver for turbulent compressible flows on unstructured meshes, with Open-FOAM implementation, Comput. Fluids 152 (2017) 14–23.

[41] V. Kolář, Vortex identification: new requirements and limitations, Int. J. Heat Fluid Flow 28 (4) (2007) 638–652.

[42] J.A. Rice, Mathematical statistics and data analysis, Cengage Learning, 2006.

[43] R.J. McCann, A convexity principle for interacting gases, Adv. Math. 128 (1) (1997) 153–179.

[44] A. Myronenko, X. Song, Point set registration: coherent point drift, IEEE Trans. Pattern Anal. Mach. Intell. 32 (12) (2010) 2262–2275.

[45] A. Ferrero, T. Taddei, L. Zhang, Registration-based model reduction of parameterized two-dimensional conservation laws, J. Comput. Phys. 457 (2022) 111068.

[46] A. Quarteroni, Numerical Models for Differential Problems, vol. 2, Springer, 2009.

[47] F.R. Menter, Two-equation eddy-viscosity turbulence models for engineering applications, AIAA J. 32 (8) (1994) 1598–1605.

[48] T.D. Economon, F. Palacios, S.R. Copeland, T.W. Lukaczyk, J.J. Alonso, SU2: an open-source suite for multiphysics simulation and design, AIAA J. 54 (3) (2016) 828–846.

[49] H. Jasak, OpenFOAM: open source CFD in research and industry, Int J. Nav. Archit. Ocean Eng. 1 (2) (2009) 89–94.

[50] MATLAB, R2022a The MathWorks Inc., Natick, Massachusetts, 2022.

[51] K. Carlberg, C. Farhat, J. Cortial, D. Amsallem, The GNAT method for nonlinear model reduction: effective implementation and application to computational fluid dynamics and turbulent flows, J. Comput. Phys. 242 (2013) 623–647.

[52] C. Farhat, T. Chapman, P. Avery, Structure-preserving, stability, and accuracy properties of the energy-conserving sampling and weighting method for the hyper reduction of nonlinear finite element dynamic models, Int. J. Numer. Methods Eng. 102 (5) (2015) 1077–1110.

[53] M. Yano, A.T. Patera, An LP empirical quadrature procedure for reduced basis treatment of parametrized nonlinear PDEs, Comput. Methods Appl. Mech. Eng. 344 (2019) 1104–1123.

[54] T. Taddei, Compositional maps for registration in complex geometries, arXiv preprint, arXiv:2308.15307, 2023.

[55] M.D. McKay, R.J. Beckman, W.J. Conover, A comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 42 (1) (2000) 55–61.

[56] E. Luke, E. Collins, E. Blades, A fast mesh deformation method using explicit interpolation, J. Comput. Phys. 231 (2) (2012) 586–601.