# Registration-based model reduction with local reduced order bases

Ishak Tifouti, Nicolas Barral, Tommaso Taddei

► To cite this version:

Ishak Tifouti, Nicolas Barral, Tommaso Taddei. Registration-based model reduction with local reduced order bases. AIAA Scitech 2025 -, AIAA, Jan 2025, Orlando (FL), United States. hal-05065511

HAL Id: hal-05065511

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

Submitted on 13 May 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
# Registration-based model reduction with local reduced order bases

Ishak Tifouti\*, Nicolas Barral\d, Tommaso Taddei\dd

## Abstract

We propose a localized automated nonlinear model reduction framework for rapid and reliable solution of parameterized shock-dominated flows. Our formulation exploits the adaptive procedure of [Barral et al, JCP, 2024] and a clustering technique to devise a piecewise Lagrangian approximation of the solution to the parametric problem: the application of clustering is designed to cope with parameter-induced shock-topology changes that hinder the effectiveness of standard (monolithic) Lagrangian approximations. We rely on (i) metric-based mesh adaptation to generate an accurate mesh for a range of parameters, (ii) parametric registration for the computation of a bijection $\Phi$ that tracks moving features of the solution field, and (iii) projection-based model reduction to rapidly and reliably estimate the mapped solution; finally, we develop (iv) a clustering technique to partition the parameter domain in subregions where the shock topology does not change. We present numerical investigations for a two-dimensional supersonic inviscid flow past a Gaussian bump to illustrate the many features of the methodology: in more detail, we display the performance of three clustering techniques, and we compare our approach with monolithic Lagrangian approximations.

## 1 Introduction

The growing demand for rapid and reliable solvers for real-time and many-query problems in computational fluid dynamics (CFD) has spurred interest in the development of model order reduction (MOR) methods for parametric shock-dominated flows. However, major advances are still required to generate effective reduced-order models (ROMs) for this class of problems, due to the inadequacy of linear approximation methods and the need for parameter-dependent adaptive meshes. During the past decade, Lagrangian (or registration-based) MOR methods have shown their effectiveness in approximating shock-dominated flows, thanks to their ability to track coherent features of the flow [20, 25, 26]. The goal of this work is to combine

<footnote>\*Inria team CARDAMOM, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, France; ishak.tifouti@inria.fr.</footnote>

<footnote>†Inria team CARDAMOM, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, France; nicolas.barral@inria.fr</footnote>

<footnote>‡Inria team MEMPHIS; Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251, France. tommaso.taddei@inria.fr.</footnote>
Lagrangian methods with parameter-domain partitioning (clustering) to enable the approximation of shock-dominated flows with topology changes.

We denote by $\mu$ the vector of model parameters in the region $\mathcal{P} \subset R^p$ and by $\Omega \subset R^d$ the $d$-dimensional computational domain; given the parametric field $w: \Omega \times \mathcal{P} \to R$, we note $w_{\mu} = w(\cdot; \mu): \Omega \to R$. We consider parametric steady conservation laws of the form

$$ \nabla \cdot F_{\mu}(q_{\mu}^{\text{true}}) = S_{\mu}(q_{\mu}^{\text{true}}) \quad \text{in } \Omega, \qquad (1) $$

completed with suitable boundary conditions; here, $q_{\mu}^{\text{true}}: \Omega \to R^D$ is the vector of $D$ state variables, $F: R^D \times \mathcal{P} \to R^{D \times d}$ is the physical flux and $S: R^D \times \mathcal{P} \to R^D$ is the source term. We denote by $\mathcal{M} = \{q_{\mu}^{\text{true}}: \mu \in \mathcal{P}\}$ the solution set, which collects the solution to (1) for all parameters in the prescribed parameter range.

Model reduction methods seek approximations in a low-dimensional linear or nonlinear space that is specifically constructed for the parametric problem at hand. MOR techniques aim to approximate the solution to (1) as $\hat{q}_{\mu} = Z_{\mu}(\hat{\alpha}_{\mu})$ where $Z: R^n \times \mathcal{P} \to \mathcal{X} := L^2(\Omega; R^D)$ is the reduced-order basis (ROB) that is built during an expensive offline (or training) stage, while $\hat{\alpha}: \mathcal{P} \to R^n$ is the vector of generalized coordinates (or latent variables) that are computed during the online (or prediction) stage for any new parameter value $\mu \in \mathcal{P}$, by solving a suitable ROM. Linear MOR methods rely on a linear parameter-independent ROB $Z$. The Kolmogorov $n$-width measures the worst-case error for the best rank-$n$ linear ROB $Z$ over the solution set $\mathcal{M}$ [28]. Projection-based MOR (PMOR) refers to a class of techniques that resort to the projection of the high-fidelity (HF) equations onto a low-dimensional test space to define the ROM for the generalized coordinates.

The past decade has witnessed a surge in the development of nonlinear approximation methods for model order reduction. These techniques are designed to overcome — or at least mitigate — the well-known issue of the slow decay of the Kolmogorov $n$-width for advection-dominated problems [27]. To provide concrete references, we recall the piecewise (or localized) methods [4, 17], quadratic and ANN-augmented methods [7, 8], the previously-mentioned Lagrangian methods, shifted (or multi-frame) methods [31, 30], and fully-nonlinear methods based on autoencoders [24, 23]. Below we focus on localized methods and Lagrangian methods.

Localized (or dictionary-based [6]) methods generate multiple local approximation bases that are dedicated to different subsets of the solution set. Given the partition $\{\mathcal{P}_i\}_{i=1}^k$ of $\mathcal{P}$, we introduce the dictionary of ROBs $\{Z_1, \dots, Z_k\}$ and we define the nonlinear ansatz:

$$ Z_{\mu}(\alpha) := \sum_{i=1}^{k} Z_{i}(\alpha) \chi_{\mathcal{P}_{i}}(\mu), \qquad (2) $$

where $\chi_{\mathcal{P}_i}(\mu) = 1$ if $\mu \in \mathcal{P}_i$ and zero otherwise. The partition $\{\mathcal{P}_i\}_{i=1}^k$ can be built in a purely data-driven fashion using clustering techniques (e.g., [4]) or by combining monolithic MOR techniques driven by a posteriori error estimators with data-driven (hierarchical) techniques [17]: in the former class of methods, clustering is applied to the dataset of simulations — i.e., in the state space — to ensure adherence to the physical problem of interest. Early works on localized approaches resort to the $k$-means method based on the Euclidean distance [4]; more recent efforts consider tree-based approaches [3, 18], or the $k$-medoids algorithm based on a more sophisticated dissimilarity measure [13]. In particular, for aerodynamic applications, Dupuis and collaborators [16] considered a clustering procedure, referred to as local decomposition method, that exploits a shock sensor and a Gaussian mixture model (GMM) to separate the solutions with a subsonic behavior from the solutions associated to the transonic regime.
Lagrangian approximations rely on a coordinate transformation to track specific features of the solution field, and ultimately simplify the task of model reduction. They consider the ansatz

$$Z_{\mu}(\alpha) = \sum_{i=1}^{n} \alpha_{i} \zeta_{i} \circ \Phi_{\mu}^{-1}, \qquad (3)$$

where $\alpha = \text{vec}(\alpha_1, \dots, \alpha_n)$, $\zeta_1, \dots, \zeta_n \in \mathcal{X}$ are linearly-independent fields, and $\Phi_{\mu}: \Omega \to \Omega$ is a parameter-dependent bijection (one-to-one and onto) in $\Omega$ that depends on an additional set of generalized coordinates, $\Phi_{\mu} := \text{N}(\hat{a}_{\mu})$ with $\text{N}: R^m \to \text{Lip}(\Omega; R^d)$. The problem of finding the parametric bijection $\Phi: \Omega \times \mathcal{P} \to \Omega$ is referred to as registration problem [32]. PMOR methods based on Lagrangian approximations exploit projection to estimate the state for out-of-sample parameters: in [12, 29, 32], the authors rely on a regression method to estimate the mapping coefficients $\hat{a}_{\mu}$ and on projection to estimate the solution coefficients $\hat{\alpha}_{\mu}$; instead, Mirhoseini and Zahr [25] rely on residual minimization for both mapping and solution coefficients.

The interest in Lagrangian methods stems from the observation that for relevant problems in continuum mechanics coherent structures that are troublesome for linear approximations such as shocks and shear layers vary smoothly with the parameter. The action $\Phi_{\mu}^{-1}$ should hence track sharp features (e.g., shocks or shear layers) of the solution field and ultimately enable the approximation of the mapped field $q_{\mu}^{\text{true}} \circ \Phi_{\mu}$ using a low-rank linear approximation $\alpha \mapsto \sum_{i=1}^{n} \alpha_i \zeta_i$ and a parameter-independent mesh; note that (see, e.g., [9, 25]) the composition with $\Phi_{\mu}^{-1}$ corresponds to a mesh deformation that can be naturally handled in the finite element (FE) or finite volume (FV) framework. Lagrangian methods have been extensively applied to both inviscid [9, 25] and viscous [29] shock-dominated flows; however, the need to ensure the bijectivity of $\Phi_{\mu}$ prevents approximations of the form (3) from effectively representing parameter-induced shock-topology changes.

We propose a piecewise-Lagrangian model reduction procedure for parametric steady shock-dominated flows,

$$Z_{\mu}(\alpha) = \sum_{i=1}^{k} \left( \sum_{j=1}^{n_i} \alpha_j \zeta_{i,j} \circ \Phi_{i,\mu}^{-1} \right) \chi_{P_i}(\mu), \qquad (4)$$

where $\zeta_{i,1}, \dots, \zeta_{i,n_i} \in \mathcal{X}$ and $\Phi_{i,\mu} := \text{N}(\hat{a}_{i,\mu})$, $i = 1, \dots, k$. Our procedure exploits the adaptive MOR and mesh adaptation (MA) framework developed in [9] and a data-driven clustering technique to handle shock-topology changes. Following [9], in each element of the partition, we first update the mesh used for HF computations, then we compute the parametric bijection $\Phi$, and finally we update the ROB $\zeta_1, \dots, \zeta_n$ and the ROM; furthermore, we rely on a data-driven clustering technique to partition regions of the domain where the low-rank representation that is needed to meet a user-defined tolerance is too large. We consider two distinct clustering techniques that are new in the framework of MOR: the first method is based on $k$-SVD [2, 10], while the second method combines the $k$-medoids algorithm with a dissimilarity measure that evaluates the distance (in the sense of Wasserstein) between suitable point clouds associated with the shocks. In the numerical experiments, we further compare our clustering procedures with the more standard $k$-means method.

Mesh adaptation is applied within each element of the partition and leads to the definition of $k$ independent meshes. The combination of registration techniques for coordinate transformation — in effect, $r$-adaptation — and clustering procedures facilitates the task of MA and ultimately enables accurate approximations with meshes of moderate size. We further remark that (4) should be coupled with a classification algorithm to determine the appropriate Lagrangian approximation that should be employed for given out-of-sample parameters: in this work, we rely on a simple nearest-neighbor classifier.
The inadequacy of Lagrangian approximations to adequately represent parametric systems with multiple transport phenomena is the primary motivation for the shifted methods in [31, 30]: our piecewise approach does not require communications between different local ROBs and is hence nearly as simple to implement as Lagrangian methods; however, due to its discontinuous nature with respect to the parameter, it might be less suited to deal with unsteady problems. A thorough comparison between the two approaches is beyond the scope of the present work.

The outline of the paper is as follows. In section 2, we introduce our procedure and we clarify inputs and outputs of the five pieces — snapshot generation, mesh adaptation, registration, model reduction, and clustering — of the formulation; in section 3, we discuss the two clustering procedures; in section 4 we introduce the model problem considered in the numerical experiments; then, in section 5, we present the numerical results. Section 6 concludes the paper.

## 2 Formulation

We rely on a nodal-based discontinuous Galerkin (DG) FE discretization. We denote by $\mathcal{T}_{hf} = (\{x_j^{hf}\}_{j=1}^{N_{nd}}, \mathbf{T})$ a mesh of the domain $\Omega$ with $N_{nd}$ nodes $\{x_j^{hf}\}_j$ and connectivity matrix $\mathbf{T}$; given the bijection $\Phi: \Omega \to \mathbb{R}^d$, we use notation $\Phi(\mathcal{T}_{hf})$ to refer to the mesh with deformed nodes $\{\Phi(x_j^{hf})\}_j$ and the same connectivity $\mathbf{T}$ as $\mathcal{T}_{hf}$. We denote by $\mathcal{X}_{hf}$ the FE space associated with the mesh $\mathcal{T}_{hf}$; any FE field $u \in \mathcal{X}_{hf}$ is uniquely characterized by the pair $(\mathcal{T}_{hf}, \mathbf{u})$ where $\mathbf{u} \in \mathbb{R}^{N_{hf}}$ is the corresponding FE vector. It is easy to verify that if $(\Phi(\mathcal{T}_{hf}), \mathbf{u})$ interpolates the field $u$ in the nodes of $\Phi(\mathcal{T}_{hf})$, then $(\mathcal{T}_{hf}, \mathbf{u})$ interpolates $u \circ \Phi$ in the nodes of $\mathcal{T}_{hf}$; the latter implies that the ansatz (3) can be stored as

$$ \mu \in \mathcal{P} \mapsto (\Phi_{\mu}(\mathcal{T}_{hf}), \mathbf{Z}\hat{\alpha}_{\mu}), \qquad (5) $$

where $\mathbf{Z} \in \mathbb{R}^{N_{hf} \times n}$ is the rank-$n$ matrix associated to the ROB $\zeta_1, \dots, \zeta_n$, and $\Phi: \Omega \times \mathcal{P} \to \Omega$ is a parametric bijection.

Exploiting (5), we can store (4) as the collection of $k$ tuples that comprise (i) a mesh $\mathcal{T}_{hf}^{(i)}$, (ii) a ROB $\mathbf{Z}^{(i)} = [\zeta_{i,1}, \dots, \zeta_{i,n_i}]$ defined over $\mathcal{T}_{hf}^{(i)}$, (iii) a parametric bijection $\mu \mapsto \hat{\Phi}_{\mu}^{(i)}$, and (iv) a ROM for the generalized coordinates $\mu \mapsto \hat{\alpha}_{\mu}^{(i)} \in \mathbb{R}^{n_i}$; we further require (v) a classifier $\hat{f}^c: \mathcal{P} \to \{1, \dots, k\}$ that identifies the appropriate local approximation for any out-of-sample configuration. In formulas, we can express (4) as

$$ \mu \in \mathcal{P} \mapsto (\Phi_{\mu}^{(i)}(\mathcal{T}_{hf}^{(i)}), \mathbf{Z}^{(i)}\hat{\alpha}_{\mu}^{(i)}), \quad \text{where } i = \hat{f}^c(\mu). \qquad (6) $$

Notice that (6) implicitly defines the partition $\{\mathcal{P}_i\}_{i=1}^k$ as $\mathcal{P}_i := \{\mu \in \mathcal{P} : \hat{f}^c(\mu) = i\}$ for $i = 1, \dots, k$.

### 2.1 Adaptive procedure

Algorithm 1 extends the adaptive procedure of [9] to piecewise approximations. To ease the presentation, we first review the core elements of the methodology; below, $\mathcal{P}_{train} = \{\mu^l\}_{l=1}^{n_{train}}$ denotes a discretization of the parameter domain, while the symbol $\sim$ refers to the mapped fields, that is $\tilde{q}_{\mu} := q_{\mu} \circ \Phi_{\mu}$. Note that (i)-(v) below are applied to each cluster separately (they hence do not involve any communication among different clusters).

(i) Snapshot generation: $(\mathcal{T}_{hf}, \Phi, \mathcal{P}_{train}) \to \{q_{\mu}^{hf} : \mu \in \mathcal{P}_{train}\}$. We generate snapshots of the solution field for all values of the parameter $\mu$ in $\mathcal{P}_{train}$ based on the parametric mesh $\mu \mapsto \Phi_{\mu}(\mathcal{T}_{hf})$.
(ii) **Mesh adaptation**: {$\tilde{q}_{\mu}^{hf} : \mu \in \mathcal{P}_{\text{train}}$} $\rightarrow \mathcal{T}_{\text{hf}}$. We build a parsimonious mesh that ensures the adequate approximations of all the elements of the (mapped) snapshot set $\tilde{\mathcal{M}}_{\text{train}} := \{\tilde{q}_{\mu}^{hf} : \mu \in \mathcal{P}_{\text{train}}\}$.

(iii) **Registration**: ($\{q_{\mu}^{hf} : \mu \in \mathcal{P}_{\text{train}}\}$, $\mathcal{T}_{\text{hf}}$) $\rightarrow \Phi$. We exploit the available set of snapshots to find a parametric bijection $\Phi$ that tracks coherent, parameter-dependent structures of the solution field. Registration should ensure that the deformed mesh $\Phi_{\mu}(\mathcal{T}_{\text{hf}})$ is a proper mesh of $\Omega$ for all parameters $\mu$ in a neighborhood of the training set $\mathcal{P}_{\text{train}}$.

(iv) **Linear-subspace model reduction**: $(\mathcal{T}_{\text{hf}}, \Phi, \mathcal{P}_{\text{train}}) \rightarrow \{\mathbf{Z}, \text{ROM}\}$. We build a linear-subspace projection-based ROM for the HF mapped solution field $\mu \mapsto q_{\mu}^{\text{hf}} \circ \Phi_{\mu}$ defined over the mesh $\mathcal{T}_{\text{hf}}$.

(v) **Clustering**:<fnref n="1" /> ($\{\tilde{q}_{\mu}^{\text{hf}} : \mu \in \mathcal{P}_{\text{train}}\}$, $k$) $\rightarrow \hat{I}_{\text{train}}^{\text{c}}$. Given the snapshot set $\tilde{\mathcal{M}}_{\text{train}}$ associated with $\mathcal{P}_{\text{train}} = \{\mu^{\ell}\}_{\ell=1}^{n_{\text{train}}}$ and the integer $k$, the clustering algorithm returns the set of indices $\hat{I}_{\text{train}}^{\text{c}} = \{i_{\ell}^{c}\}_{\ell=1}^{n_{\text{train}}} \subset \{1, \dots, k\}$ such that $i_{\ell}^{c}$ is the label of the cluster to which $\tilde{q}_{\mu^{\ell}}^{\text{hf}}$ belongs.

(vi) **Classification**: {$\mu^{\ell}, i_{\ell}^{c}$}$_{\ell=1}^{n_{\text{train}}} \rightarrow \hat{f}^{\text{c}}$. Given the labeled set of parameters $\{\mu^{\ell}, i_{\ell}^{c}\}_{\ell=1}^{n_{\text{train}}}$, we apply classification to determine a partition of the entire parameter domain.

The core elements (i)-(vi) are combined in an iterative algorithm that adaptively builds and improves the local reduced order approximations, the local meshes, and the local reduced-order models.

## 2.2 Discussion

The first four building blocks are the same as in [9]; in particular, we recall that snapshot generation is accelerated by resorting to the available ROM constructed during the previous iteration of the outer loop (cf. Lines 4-17 of Algorithm 1). We refer to [9] for a detailed presentation of the algorithms for the tasks (i)-(iv). We further observe that the snapshots generated at Line 15 of Algorithm 1 can be reused at Line 6 of the subsequent iteration. We notice that the criterion employed to update the clusters (cf. Line 14) ensures that the local ROBs cannot exceed the user-defined threshold $n_{\text{max}}$: a similar criterion was proposed in [17, Algorithm 3].

Notice that the number of clusters might change during the iterative loop: for instance, if the criterion in Line 14 is always met, we find that the number of clusters at the *i*-th iteration $k^{(i)}$ satisfies $k^{(i)} = 2k^{(i-1)}$ for $i = 1, \dots, N_{\text{it}}$. Notice also that clustering and classification are hierarchical: if the parameters $\mu$ and $\mu'$ belong to two separate clusters at iteration $k$, they are indeed guaranteed to belong to different clusters for all the subsequent iterations, provided that the classification algorithm is “exact” in the training set. This decision is motivated by the observation that “moving” one datapoint from one cluster to another requires communication between different meshes and/or it requires to extrapolate the local Lagrangian approximation outside its domain of definition. In the numerical experiments, we resort to nearest-neighbor classification.

<footnote n="1">¹As discussed in the next section, the clustering procedure might either require the snapshots in the reference configuration (cf. section 3.1) or the snapshots in the physical configuration (cf. section 3.2). To ease the presentation, we do not distinguish the two cases in Algorithm 1.</footnote>
**Algorithm 1**: local adaptive training (offline) procedure.

**Initialization:**

define the mesh $\mathcal{T}_{\text{hf}}^{(0)}$,
$\Phi^{(0)} = \text{id}$ identity map,
$\mathcal{P}_{\text{train}} = \{\mu^{\ell}\}_{\ell=1}^{n_{\text{train}}}$ initial training set,
and the initial number of clusters $k^{(0)}$.

1: Snapshot generation
   $(\mathcal{T}_{\text{hf}}^{(0)}, \Phi^{(0)}, \mathcal{P}_{\text{train}}) \rightarrow \{q_{\mu}^{\text{hf}, (0)} : \mu \in \mathcal{P}_{\text{train}}\}$

2: Apply clustering and classification to define the partition
   $(\{q_{\mu}^{\text{hf}, (0)} : \mu \in \mathcal{P}_{\text{train}}\}, k^{(0)}) \rightarrow \hat{I}_{\text{train}}^{c}$
   $\{\mu^{\ell}, i_{\ell}^{c}\}_{\ell=1}^{n_{\text{train}}} \rightarrow \hat{f}^{c, (0)}$

3: **for** $i = 1, \dots, N_{\text{it}}$ **do**

4: &nbsp;&nbsp;&nbsp;&nbsp;**for** $j = 1, \dots, k^{(i-1)}$ **do**

5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $\mathcal{P}_{\text{train}}^{(j)} := \{\mu \in \mathcal{P}_{\text{train}} : \hat{f}^{c, (i-1)}(\mu) = j\}$

6: Snapshot generation
   $(\mathcal{T}_{\text{hf}}^{(i-1, j)}, \Phi^{(i-1, j)}, \mathcal{P}_{\text{train}}^{(j)}) \rightarrow \{q_{\mu}^{\text{hf}, (i, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}\}$

7: **if** $i > 1$ **then**

8: &nbsp;&nbsp;&nbsp;&nbsp;Mesh adaptation
   $\{q_{\mu}^{\text{hf}, (i, j)} \circ \Phi_{\mu}^{(i-1, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}\} \rightarrow \mathcal{T}_{\text{hf}}^{(i, j)}$

9: **else**

10: Set $\mathcal{T}_{\text{hf}}^{(i,j)} = \mathcal{T}_{\text{hf}}^{(0)}$.

11: **end if**

12: Registration
   $\{(q_{\mu}^{\text{hf}, (i, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}, \mathcal{T}_{\text{hf}}^{(i, j)}) \rightarrow \Phi^{(i, j)}\}$

13: Linear-subspace model reduction
   $(\mathcal{T}_{\text{hf}}^{(i, j)}, \Phi^{(i, j)}, \mathcal{P}_{\text{train}}^{(j)}) \rightarrow (\mathbf{Z}^{(i, j)}, \text{ROM}^{(i, j)})$

14: **if** $n_j > n_{\text{max}}$ **then**

15: &nbsp;&nbsp;&nbsp;&nbsp;Snapshot generation
   $\{q_{\mu}^{\text{hf}, (i+1, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}\}$
   $(\{\tilde{q}_{\mu}^{\text{hf}, (i+1, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}, k=2\} \rightarrow \hat{I}_{\text{train}}^{c, (i, j)}$

16: Clustering
   $\{(\tilde{q}_{\mu}^{\text{hf}, (i+1, j)} : \mu \in \mathcal{P}_{\text{train}}^{(j)}, k=2) \rightarrow \hat{I}_{\text{train}}^{c, (i, j)}\}$

17: **end if**

18: **end for**

19: Exploit $\{\mathcal{P}_{\text{train}}^{(j)}, \hat{I}_{\text{train}}^{c, (i, j)}\}_{j=1}^{k^{(i-1)}}$ to define the labeled set $\{\mu^{\ell}, i_{\ell}^{c, (i)}\}_{\ell=1}^{n_{\text{train}}}$
   and apply classification
   $\rightarrow \hat{f}^{c, (i)}$

20: **end for**

We empirically observed that the algorithm might lead to poor performance if $n_{\text{train}}/k$ is too small.
Since the number of clusters $k$ increases during the iterative loop, we propose to add training points to $\mathcal{P}_{\text{train}}$: in the numerical experiments, we add random samples in each cluster at the end of each outer-loop iteration $i$ to ensure that the cardinality of $\mathcal{P}_{\text{train}}^{(i)}$ is above a user-defined threshold $n_{\text{train}}^{\min}$. Given the classifier $\hat{f}^{c,(i)}$ and the label $j \in \{1, \dots, k^{(i)}\}$, we resort to rejection sampling to generate independent identically distributed (i.i.d.) samples for each cluster $j$:

<table><thead><tr><th colspan="2">Algorithm 2 Rejection sampling algorithm</th></tr></thead><tbody><tr><td>while true do<br/>Sample &mu; ~ Uniform(<i>P</i>).</td><td></td></tr><tr><td>if $\hat{f}^{c,(i)}(\mu) = j$ then, Return &mu;</td><td></td></tr><tr><td>end if</td><td></td></tr><tr><td>end while</td><td></td></tr></tbody></table>

Finally, we observe that our approximation is based on a non-overlapping partition of the parameter domain both for training — that is, $\mathcal{P}_{\text{train}}^{(i)} \cap \mathcal{P}_{\text{train}}^{(j)} = \emptyset$ if $i \neq j$ (cf. Line 6, Algorithm 1) — and for prediction — that is, given $\mu \in \mathcal{P}$, we only evaluate the ROM associated with the label returned by the classifier (cf. (6)). The authors of [4] advocate the use of an overlapping partition (cf. [4, Algorithm 3]) to deal with transitions from one cluster to another: while this issue is crucial for unsteady PDEs, it does not appear relevant for steady-state problems.

# 3 Clustering

In this work, we consider two clustering methods: $k$-SVD and $k$-medoids. Both techniques ensure that each element of the dataset is assigned to one cluster (*hard clustering*) and can be interpreted as generalizations of $k$-means clustering [19]. We review the key features of the methods in sections 3.1 and 3.2. Both methods are applied in the state space — i.e., they are applied to the solution snapshots — to ensure the adherence to the physical problem of interest; however, as discussed in section 3.3, we postprocess the results with a linear classification algorithm in the parameter space to ensure that the clusters are connected and have a simple shape — the latter property is important for the classification step performed in the prediction phase (cf. (6)) and also for parametric registration, which relies on regression for generalization (see, e.g., [9]).

## 3.1 $k$-SVD

Given the integers $n$ and $k$, the $k$-SVD algorithm [2] aims to minimize the sum of the average projection errors associated with the SVD linear space of dimension $n$ in each cluster. We denote by $\{\tilde{q}_\ell\}_{\ell=1}^{n_{\text{train}}} \subset \mathcal{X}_{\text{hf}}$ the training set of simulations; we introduce the Hilbert norm $\| \cdot \|$ induced by the inner product $(\cdot, \cdot)$ in $\mathcal{X}_{\text{hf}}$, $\| \cdot \| = \sqrt{(\cdot, \cdot)}$; given the linear space $\mathcal{Z} \subset \mathcal{X}_{\text{hf}}$, we define the orthogonal projection operator $\Pi_{\mathcal{Z}} : \mathcal{X}_{\text{hf}} \to \mathcal{Z}$. We further recall that, given $n < n_{\text{train}}$, proper orthogonal decomposition<fnref n="2" /> (POD, [33])

<footnote n="2">²In the finite-dimensional setting, POD is a weighted SVD. In the original paper [2], the authors used the Euclidean norm, while here we prefer to exploit the $L^2(\Omega)$ norm.</footnote>
identifies the space $\mathcal{Z}_n^{\text{pod}}$ of dimension $n$ that minimizes the total projection error

$$ \varepsilon_n (\{\tilde{q}_\ell\}_{\ell=1}^{n_{\text{train}}}) := \min_{\mathcal{Z} \subset \mathcal{X}_{\text{hf}}, \dim(\mathcal{Z})=n} \sum_{\ell=1}^{n_{\text{train}}} \|\tilde{q}_\ell - \Pi_{\mathcal{Z}} \tilde{q}_\ell\|^2 = \sum_{\ell=1}^{n_{\text{train}}} \|\tilde{q}_\ell - \Pi_{\mathcal{Z}_n^{\text{pod}}} \tilde{q}_\ell\|^2. $$

Exploiting the previous definitions, k-SVD seeks a vector of labels $\mathbf{I} = \text{vec}(\mathbf{I}_1, \dots, \mathbf{I}_{n_{\text{train}}})$ where $\mathbf{I}_\ell \in \{1, \dots, k\}$ that minimizes

$$ f(\mathbf{I}) := \sum_{i=1}^{k} \varepsilon_n (\{\tilde{q}_\ell : \mathbf{I}_\ell = i\}). \qquad (7) $$

Similarly to the k-means algorithm, k-SVD alternates between an assignment step in which each datapoint is assigned to the cluster with the lowest projection error, and an update step in which we recalculate the SVD spaces for each cluster. As discussed in [10], the iterative procedure is guaranteed to converge to a local minimum. We provide the Algorithm in Appendix A for completeness.

The performance of k-SVD is highly sensitive to the choice of the initial guess of the alternating procedure; in the numerical experiments, we use the k-means algorithm with the initialization described in [5] to generate the initial guess. k-SVD is inadequate to capture shock-topology changes since it is based on linear dimensionality reduction: for this reason we conjecture that its application at Line 2 of Algorithm 1 might lead to a poor partitioning of the parameter space. However, in the subsequent iterations, the application of k-SVD to the mapped snapshots is effective to separate the effect of transport of coherent structure (for which registration is effective) from the effects associated with topology changes (for which registration is ineffective).

## 3.2 k-medoids

In contrast to k-means, the k-medoids algorithm [22] chooses data points as centers of each cluster: this feature of the method typically produces more interpretable results. Furthermore, the algorithm takes as input a dissimilarity matrix $\mathbf{D} \in \mathbb{R}^{n_{\text{train}} \times n_{\text{train}}}$ such that $(\mathbf{D})_{i,j}$ is the dissimilarity between the $i$-th and the $j$-th datapoint: it can hence cope with arbitrary dissimilarity measures. In this work, we exploit the latter feature to define a dissimilarity measure based on a Gaussian mixture model (GMM) of the set of the discontinuities of the solution field. We anticipate that the dissimilarity measure considered is a pseudo-metric: it satisfies non-negativity, symmetry, and the triangle inequality (but may allow zero distance between distinct points). We further expect that such a pseudo-metric still enables k-medoids to group points in a way that minimizes intra-cluster distances.

Given the state $q$ of conserved variables, we denote by $u$ the velocity field, by $p$ the pressure, by $a$ the sound speed, and by $Ma$ the Mach number. Given the mesh $\mathcal{T}_{\text{hf}}$, we define the elements $\{\mathbf{D}_\ell\}_{\ell=1}^{N_e}$ with centroids $\{x_\ell^{\text{hf,c}}\}_{\ell=1}^{N_e}$. Given the scalar or vector-valued function $f$, we define the element averages $\tilde{f}_\ell := \frac{1}{|\mathbf{D}_\ell|} \int_{\mathbf{D}_\ell} f(x) dx$ for $\ell = 1, \dots, N_e$; if $f$ is scalar, we also define the maximum $\tilde{f}_\ell^{\text{max}} := \max_{x \in \mathbf{D}_\ell} f(x)$ for $\ell = 1, \dots, N_e$. Then, we introduce the shock sensor $\eta[q] \in \mathbb{R}^{N_e}$ such that

$$ (\eta[q])_\ell = \begin{cases} \frac{1}{|\mathbf{D}_\ell|} \int_{\mathbf{D}_\ell} \frac{(-\nabla \cdot u)_+}{\sqrt{(\nabla \cdot u)^2 + |\nabla \times u|^2 + a^2}} dx & \text{if } \overline{\text{Ma}}_\ell^{\text{max}} > 1, \overline{\nabla p}_\ell \cdot \overline{u}_\ell > 0 \\ 0 & \text{otherwise} \end{cases} \quad \ell = 1, \dots, N_e. \quad (8) $$

The first term in (8) is a variant of the Ducros sensor [15], whose expression is justified by the fact that shocks occur in strong compression regions where rotation effects are negligible compared to compression
regions; the condition $Ma_{\ell}^{\max} > 1$ ensures that the sensor vanishes in the subsonic region of the domain; finally, the condition $\nabla p_{\ell} \cdot \bar{u}_{\ell} > 0$ is motivated by the observation that pressure increases across a shock wave. Finally, we define the point cloud associated with the state

$$C[q] = \{x_{\ell}^{\text{hf,c}} : (\eta[q])_{\ell} \ge \text{tol}\}, \quad (9)$$

where tol is a user-defined threshold that is set equal to 0.7 in the numerical experiments. The set $C[q]$ describes the discontinuity set of the function q.

Following [21], we interpret the elements of $C[q]$ as i.i.d. samples from an unknown distribution $\rho[q]$; we approximate $\rho[q]$ using a GMM. We recall that a GMM is a random variable with probability density function

$$\rho_{\text{gmm}}(x; \theta) = \sum_{j=1}^{M} \pi_j \rho_{\text{gm}}(x; m_j, \Sigma_j), \quad (10)$$

where $\rho_{\text{gm}}(\cdot; m, \Sigma)$ is the probability density function associated with a multivariate Gaussian random variable with mean $m$ and covariance $\Sigma$, and $0 \le \pi_1, \dots, \pi_M \ge 1$ are the mixture coefficients such that $\sum_j \pi_j = 1$. Given the point cloud $\{x_k^{\text{hf,c}} : (\eta[q])_k \ge \text{tol}\}$, we apply the expectation-maximization algorithm reviewed in [19, Chapter 8.5] to build the GMM; in the numerical experiments, we exploit the Matlab function fitgmdist and we select the number mixture modes $M$ in (10) using the Aikake information criterion (see, e.g., [19, Chapter 7.5]).

Next, we exploit the densities $\rho_{\text{gmm}}^{\ell} := \rho_{\text{gmm}}[q_{\mu^{\ell}}^{\text{hf}}]$ for $\ell = 1, \dots, n_{\text{train}}$ to compute the dissimilarity matrix $\mathbf{D}$. The Wasserstein distance measures the cost of “transporting” a given probability distribution into another: it is hence well-suited to measure the distance between different shock topologies. Since the exact computation of the Wasserstein distance is expensive, we propose to use the approximate Wasserstein distance that is valid for GMM and was proposed in [11, 14]. First, we recall the explicit formula for the 2-Wasserstein distance between two Gaussian random variables with mean and variance $(m_0, \Sigma_0)$ and $(m_1, \Sigma_1)$:

$$W_2^2(\rho_{\text{gm}}^0, \rho_{\text{gm}}^1) = \sqrt{\|m_0 - m_1\|_2^2 + \text{trace} \left(\Sigma_0 + \Sigma_1 - 2 \left(\Sigma_0^{1/2} \Sigma_1 \Sigma_0^{1/2}\right)^{1/2}\right)}. \quad (11)$$

Then, given $\rho_{\text{gmm}}^0 = \sum_{i=1}^{M_0} \pi_i^0 \rho_{\text{gm}}^{0,i}$ and $\rho_{\text{gmm}}^1 = \sum_{j=1}^{M_1} \pi_j^1 \rho_{\text{gm}}^{1,j}$, we define the approximate Wasserstein distance of [11, 14] as

$$MW_2(\rho_{\text{gmm}}^0, \rho_{\text{gmm}}^1) = \min_{\mathbf{w} \in \mathcal{R}^{M_0 \times M_1}} \sum_{i=1}^{M_0} \sum_{j=1}^{M_1} (\mathbf{w})_{i,j} W_{i,j} \quad \text{s.t.} \quad \left\{ \begin{array}{ll} \sum_{j=1}^{M_1} (\mathbf{w})_{i,j} = \pi_i^0 & i = 1, \dots, M_0 \\ \sum_{i=1}^{M_0} (\mathbf{w})_{i,j} = \pi_j^1 & j = 1, \dots, M_1 \end{array} \right. \quad (12)$$

where $W_{i,j} := W_2^2(\rho_{\text{gm}}^{0,i}, \rho_{\text{gm}}^{1,j})$ for $i = 1, \dots, M_0$ and $j = 1, \dots, M_1$. Note that computation of (12) requires the solution to a linear programming problem of size $M_0 \cdot M_1$.

Algorithm 3 summarizes the whole procedure. Notice that, while in k-SVD the clustering procedure takes as input the mapped snapshots to facilitate the application of POD, the procedure takes as input the snapshots in the deformed mesh. Notice also that since the snapshots are only used to generate the point clouds, we can readily consider snapshots defined over different meshes.
**Algorithm 3 : k-medoids clustering.**

1: Compute the sensors $\eta[q_{\mu^1}^{hf}], \dots, \eta[q_{\mu^{n_{train}}}^{hf}]$ in (8) and define the point clouds $C[q_{\mu^1}^{hf}], \dots, C[q_{\mu^{n_{train}}}^{hf}]$ in (9).

2: Compute the GMM $\rho_{gmm}^k := \rho_{gmm}[q_{\mu^k}^{hf}]$, for $k = 1, \dots, n_{train}$.

3: Compute the distances $(\mathbf{D})_{i,j} = MW_2(\rho_{gmm}^i, \rho_{gmm}^j)$ in (12) for $i, j = 1, \dots, n_{train}$.

4: Apply the k-medoids algorithm [22], which takes as input $\mathbf{D}$ and returns the labels $\hat{I}_{train}^c$.

## 3.3 Regularization

In order to facilitate the task of classification, it is important to ensure that each cluster has a “simple shape” in the parameter domain. Since clustering is applied in the state space, we propose to post-process the labels by applying a linear classifier based on the labels returned by the clustering procedure.

Algorithm 4 summarizes the computational procedure: in the numerical experiments, we consider logistic regression when $k = 2$ (see the documentation of the Matlab function `fitglm`) and one-vs-one support vector machine for $k > 2$ (see the documentation of the Matlab function `fitcecoc`). We use notation $\hat{I}_{train}^{c,raw}$ to indicate the labels predicted by the clustering procedure and notation $\hat{I}_{train}^c$ to indicate the final labels that are returned by the algorithm. Figure 1 provides a pictorial explanation of the post-processing stage for $k = 2$ and a two-dimensional parameter domain.

**Algorithm 4 : regularized clustering.**

1: Apply the clustering procedure (cf. Algorithms 5 or 3) $\{\tilde{q}_{\mu}^{hf} : \mu \in \mathcal{P}_{train}\} \to \hat{I}_{train}^{c,raw}$.

2: Train a linear classifier $\{\mu^\ell, i_\ell^{c,raw}\} \to \tilde{f}_{\text{lin}}^c$.

3: Return the labels $\hat{I}_{train}^c := \{i_\ell^c = \tilde{f}_{\text{lin}}^c(\mu^\ell) : \ell = 1, \dots, n_{train}\}$.
<figure><img src="image_2.png" /><figcaption>Figure 1: illustration of the regularization process. (a) Labels returned by the clustering procedure and corresponding linear separating hyper-plane. (b) Labels after post-processing</figcaption></figure>

## 4 Model problem

We consider a two-dimensional inviscid flow past a Gaussian bump. We introduce the domain $\Omega = \{x \in (-1.5, 1.5) \times (0, 0.8) : x_2 > he^{-25x_1^2}\}$ where $h > 0$ is a given parameter. We denote by $q = \text{vec}(\rho, \rho u, E)$ the vector of conserved variables, where $\rho$ is the density, $u$ is the velocity and $E$ is the total energy; we further introduce the pressure $p = (\gamma - 1)(E - \frac{1}{2}\rho|u|^2)$, where $\gamma > 0$ is the ratio of specific heats. We consider the steady compressible Euler equations in conservative form:

$$ \nabla \cdot F(q) = 0, \quad \text{where} \quad q = \begin{bmatrix} \rho \\ \rho u \\ E \end{bmatrix}, \quad F(q) = \begin{bmatrix} \rho u^\top \\ \rho u u^\top \\ (E+p)u^\top \end{bmatrix}, \qquad (13) $$

completed with wall boundary conditions on top and bottom boundaries, supersonic inlet condition at the left boundary and supersonic outlet condition at the right boundary — the symbol $(\cdot)^\top$ denotes the transposition operator. We express the free-stream field $q_\infty$ in terms of the Mach number $Ma_\infty$,

$$ T_\infty = 1, \quad p_\infty = \frac{1}{\gamma}, \quad \rho_\infty = 1, \quad u_\infty = Ma_\infty e_1. $$

Finally, we introduce the parameter vector $\mu = [h, Ma_\infty]$ and the parameter region $P = [0.05, 0.065] \times [1.6, 3]$. Note that the computational domain $\Omega$ 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 refer to [9] for the details.

Figure 2 shows the behavior of the Mach field for four parameter values in $P$. We notice that parameter variations induce shock topology changes, which fundamentally undermine the effectiveness of Lagrangian methods.
<figure><img src="image_3.png" /><figcaption>Figure 2: model problem; behavior of the Mach number for several parameter values. (a) $\mu = [0.06, 1.69]$. (b) $\mu = [0.058, 2.02]$. (c) $\mu = [0.065, 2.41]$. (d) $\mu = [0.065, 3]$.</figcaption></figure>

# 5 Numerical results

We investigate the performance of the piecewise Lagrangian MOR method developed in this work for the model problem of section 4. We consider a P2 discontinuous Galerkin (DG) discretization of the conservation laws with dilation-based artificial viscosity and we resort to pseudo-transient continuation (PTC) to solve the nonlinear system of equations. We rely on metric-based mesh adaptation driven by the Hessian of the Mach number; furthermore, we consider metric intersection to combine metrics associated with different snapshots; the implementation is based on the software `mmg2d` [1]. We consider a least-square Petrov-Galerkin (LSPG) ROM for the generalized coordinates; the ROB and the ROM are built adaptively using the weak-greedy algorithm with termination condition $tol_{gr} = 10^{-3}$; we also prescribe the maximum size of the ROB equal to $n_{max,greedy} = 50$. In Algorithm 1, we consider an initial parameter domain $\mathcal{P}_{train}$ with $10 \times 10$ equispaced parameters; we set the initial number of clusters equal to $k^{(0)} = 3$, and we split the cluster when the ROB is of dimension larger than $n_{max} = 25$ (cf. Line 14 of Algorithm 1). We refer to [9] for a thorough description of the many elements of the formulation. Simulations are performed in Matlab R2019a based on an in-house code, and executed over a dedicated computing machine with RAM 4 GB, 2x 32-core AMD Zen2 EPYC 7452 @ 2.35 GHz.

## 5.1 Performance analysis

### 5.1.1 Performance metrics

We assess quantitative performance based on $n_{test} = 20$ out-of-sample parameters $\mathcal{P}_{test} = \{\mu_{test}^j\}_{j=1}^{n_{test}}$ with $\mu_{test}^1, \dots, \mu_{test}^{n_{test}} \stackrel{\text{iid}}{\sim} \text{Uniform}(\mathcal{P})$. A k-nearest neighbor classifier is applied to $\mathcal{P}_{test}$ to determine the appropriate local reduced model for each testing parameter. For each $\mu \in \mathcal{P}_{test}$, we report the HF $L^2$
error $E_{\mu}^{hf}$, the total enthalpy error $E_{\mu}^{\infty}$ and the sub-optimality index $\eta_{\mu}^{hf}$ such that

$$E_{\mu}^{hf} = \frac{\|q_{\mu}^{hf} - \hat{q}_{\mu}^{hf}\|_{L^2(\Omega_{\mu})}}{\|q_{\mu}^{hf}\|_{L^2(\Omega_{\mu})}}, \qquad (14)$$

$$E_{\mu}^{\infty} = \frac{\|H_{\text{tot},\mu}^{\text{true}} - \hat{H}_{\text{tot},\mu}\|_{L^2(\Omega_{\mu})}}{\|H_{\text{tot},\mu}^{\infty}\|_{L^2(\Omega_{\mu})}}, \qquad (15)$$

$$\eta_{\mu}^{hf} = \frac{\|q_{\mu}^{hf} - \hat{q}_{\mu}^{hf}\|_{L^2(\Omega_{\mu})}}{\min_{\zeta \in \mathcal{Z}_n} \|q_{\mu}^{hf} - \zeta \circ \Phi_{\mu}^{-1}\|_{L^2(\Omega_{\mu})}}, \quad \text{where } \mathcal{Z}_n := \text{span}\{\zeta_i\}_{i=1}^n. \qquad (16)$$

The relative error $E_{\mu}^{hf}$ measures the accuracy of the reduced-order estimate with respect to the HF model employed for training — it is hence a measure of the overall ability of the MOR procedure to approximate the truth model of the PDE. The total enthalpy error (15) measures the accuracy of the state estimate with respect to the exact solution to the PDE, in terms of enthalpy preservation. Finally, $\eta_{\mu}^{hf}$ measures the extent to which the LSPG projection scheme is suboptimal compared to the best-fit error.

## 5.1.2 Results

We compare performance of five different methods: (i) the global Lagrangian method of [9]; the piecewise Lagrangian method with (ii) k-SVD clustering, (iii) k-medoids clustering, (iv) k-means clustering based on the Euclidean distance, and (v) k-means clustering based on the cosine distance (“angular k-means”). We resort to the Matlab function `kmeans` to implement k-means with the two distinct distances. To facilitate the comparison between the global and the piecewise approaches, the global Lagrangian ROM is obtained by forcing $k = 1$ in Algorithm 1 and by considering the same input parameters — i.e., initial mesh, initial training set.

Figure 3 illustrates the behavior of slices of the Mach field over the bump for $x_1 \in [0.4, 0.9]$, over the bump for five parameter values in $\mathcal{P}_{\text{test}}$. The results are compared to a HF solution computed on an adapted fine mesh with approximately 19,000 elements: the HF solution for each parameter is derived from a distinct mesh that is obtained at the end of a four-step mesh adaption loop to achieve high accuracy. The comparison demonstrates a strong agreement between the reduced solutions (dotted lines) and the high-fidelity solutions (solid lines) for both local and global reduction approaches. Qualitatively, all local methods and the global approach show comparable performance.

Table 1 shows the offline cost of each approach. We provide the number of clusters, the average size of the ROBs in each cluster, the number of mesh elements, the total offline cost, and the number of PTC iterations of the Newton solver. We observe that the global Lagrangian method reaches the prescribed maximum basis size and is considerably more expensive than the other four methods. We further observe that the cost of the piecewise methods are comparable. Interestingly, k-means with Euclidean distance is the only algorithm that returns four clusters. Offline costs are dominated by the solution the HF problems, which account for roughly 50% of the total offline costs; k-means and k-SVD clustering account for approximately 0.7% – 0.8% of the offline costs (the cost includes snapshot generation at Line 15 of Algorithm 1); on the other hand, k-medoids is slightly more expensive as it accounts for 2% of the offline costs. We conclude that the cost of clustering and classification is a small fraction of the offline training stage.
<figure><img src="image_4.png" /><figcaption>Figure 3: performance analysis; slices of the Mach field for five parameter values. Dotted lines indicate the predictions of the ROM; continuous lines correspond to the HF approximation. (a) global Lagrangian. (b) *k*-SVD. (c) *k*-medoids. (d) *k*-means with Euclidean distance. (e) *k*-means with cosine distance.</figcaption></figure>

<table><caption>Table 1: performance analysis; offline (training) stage. The numbers in brackets indicate minimum and maximum value among clusters.</caption><thead><tr><th>Offline</th><th>nbr clusters</th><th>ROB size (avg)</th><th>Mesh complexity (avg)</th><th>offline costs</th><th>PTC iterations (avg)</th></tr></thead><tbody><tr><td>Lagrangian (global)</td><td>1</td><td>50</td><td>8503</td><td>02d-22h:55m</td><td>14.18</td></tr><tr><td>*k*-SVD</td><td>5</td><td>20.2 (19-22)</td><td>8119.2 (8061-8195)</td><td>01d-06h:24m</td><td>17.4 (11.57-33.59)</td></tr><tr><td>*k*-medoids</td><td>5</td><td>21.4 (18-25)</td><td>8038.6 (7936-8088)</td><td>01d-05h-03m</td><td>13.6 (11.65-17.60)</td></tr><tr><td>*k*-means</td><td>4</td><td>23.75 (21-27)</td><td>7972.25 (7994-8062)</td><td>01d-02h:22m</td><td>14.9 (11.74-19.38)</td></tr><tr><td>Angular *k*-means</td><td>5</td><td>20.4 (17-25)</td><td>8071.6 (7935-8175)</td><td>01d-03h:16m</td><td>15.0 (10.86-19.00)</td></tr></tbody></table>

Table 2 summarizes the online performance of the global and the piecewise ROMs. First, we notice that the $L^2$ error of the global ROM is roughly one-order-of-magnitude higher than the errors of the piecewise ROMs due to the fact that the greedy algorithm exits when the ROB size reaches $n_{max,greedy} = 50$. We further notice that the online prediction cost is roughly tenfold higher. On the other hand, the piecewise ROMs exhibit comparable performance, with *k*-SVD showing a slight advantage over the other methods.
<table><caption>Table 2: performance analysis; online (prediction) stage.</caption><thead><tr><th>Online</th><th>online cost (avg)</th><th>L<sup>2</sup> error (avg)</th><th>entr. error (avg)</th><th>suboptimality (avg)</th></tr></thead><tbody><tr><td>Lagrangian (global)</td><td>2.7849s</td><td>3.65e-03</td><td>2.36e-03</td><td>1.78</td></tr><tr><td>k-SVD</td><td>0.2859s</td><td>3.20e-04</td><td>2.46e-03</td><td>1.20</td></tr><tr><td>k-medoids</td><td>0.3565s</td><td>5.18e-04</td><td>2.47e-03</td><td>1.21</td></tr><tr><td>k-means</td><td>0.3679s</td><td>4.76e-04</td><td>2.46e-03</td><td>1.62</td></tr><tr><td>Angular k-means</td><td>0.3210s</td><td>4.73e-04</td><td>2.44e-03</td><td>1.38</td></tr></tbody></table>

## 5.2 Clustering analysis

Figure 4 shows the partitions of the parameter domain after four iterations of Algorithm 1. Interestingly, all four methods show a tendency to split the parameter domain to capture variations of the free-stream Mach number. Three clustering (k-SVD, k-medoids, and k-means with cosine distance) techniques show a cluster with a small area that corresponds to the lowest values of the Mach number $Ma_∞$. Furthermore, the same three clustering procedures return a piecewise approximation with five clusters, while k-means returns a piecewise approximation with four clusters.

<figure><img src="image_5.png" /><figcaption>Figure 4: The final partitions for each clustering method</figcaption></figure>

Figures 5 and 6 show representative Mach fields associated with each element of the partition. We notice that the k-means algorithm is not able to separate solution fields with different shock topologies, while k-SVD does. This is due to the fact that k-means clustering is sensitive to variations in the solution absolute norm.
<figure><img src="image_6.png" /><figcaption>Figure 5: Partition obtained with the *k*-means algorithm.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Partition obtained with the *k*-SVD algorithm.</figcaption></figure>

The *k*-medoids algorithm differs from the other clustering methods discussed in this work as it operates on a point-based approach instead of directly operating on the solution fields. Clearly, the performance of *k*-medoids strongly depends on the ability of the sensor to identify the regions of the spatial domain where the shocks occur. Figure 7 shows the points (9) selected through the vehicle of the Ducros sensor for four distinct parameter values: we notice that the sensor is able to correctly capture the shock topology.
<figure><img src="image_8.png" /></figure>

(a) $\mu = [0.06, 1.69]$

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

(b) $\mu = [0.058, 2.02]$

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

(c) $\mu = [0.065, 2.41]$

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

(d) $\mu = [0.065, 3]$

<caption>Figure 7: Point clouds (9) selected by the Ducros sensor (8) for four parameter values.</caption>

## 5.3 Updates of the partition

We analyze the evolution of the partitioning process based on the updating criterion specified at Line 14 of Algorithm 1. Figures 8-11 display the partitions obtained across the four iterations of the algorithm for the different clustering procedures. Interestingly, *k*-SVD is the only algorithm that identifies the partition with five clusters after three iterations of the adaptive procedure, which suggests its superiority to capture shock topology changes for relatively coarse meshes.

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

(a) $it = 1$ and $it = 2$

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

(b) $it = 3$ and $it = 4$

<caption>Figure 8: Evolution of the partition of the domain *P* using *k*-SVD.</caption>
<figure><img src="image_14.png" /><figcaption>Figure 9: Evolution of the partition of the domain $\mathcal{P}$ using $k$-medoids.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 10: Evolution of the partition of the domain $\mathcal{P}$ using $k$-means.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 11: Evolution of the partition of the domain $\mathcal{P}$ using $k$-means with cosine distance.</figcaption></figure>

# 6 Conclusions

We developed a piecewise Lagrangian model reduction for shock-dominated flows with parameter-induced topology changes. The approach combines the adaptive procedure of [9] with a clustering technique to cope with parameter-induced topology changes that cannot be effectively captured by a single coordinate transformation. We discussed the performance of four clustering techniques — $k$-SVD, $k$-medoids driven by a dissimilarity measure based on a shock sensor, $k$-means and $k$-means based on cosine distance.

The numerical results for a parametric inviscid flow past a Gaussian bump clearly showed the superiority of piecewise Lagrangian approximations over monolithic (global) Lagrangian approximations both in terms of offline costs and in terms of online costs. Among the clustering techniques considered, we found $k$-SVD to be extremely effective to build appropriate partitions of the parameter domain for relatively-coarse approximations. On the other hand, our approach based on $k$-medoids has the benefit of being able to cope with snapshot sets defined over different meshes.

# Acknowledgments

The authors acknowledge the support provided by Inria through the exploratory action program (project title: *Adaptive Meshes for Model Order Reduction*, AM2OR). T.T. acknowledges the support by European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Actions, grant agreement 872442 (ARIA).
A Algorithms

Algorithm 5 k-SVD: ({\tilde{q}_j : j = 1, \dots, n_{\text{train}}}, k, n) \to \hat{I}^c

Initialization:

define the initial partition $\hat{I}^{c,(0)}$
compute the POD spaces $\mathcal{Z}^{(0,i)} := \text{POD}(\{\tilde{q}_j : \hat{I}_j^{c,(0)} = i\}, n, \|\cdot\|)$ for $i=1, \dots, k$.

1: **for** $k = 1, \dots, \maxit$ **do**

2: Set $\hat{I}_i^{c,(k)} := \arg\min_{j=1,\dots,k} \|\tilde{q}_i - \Pi_{\mathcal{Z}^{(k-1,j)}} \tilde{q}_i\|$ for $i=1, \dots, n_{\text{train}}$.

3: Compute the POD spaces $\mathcal{Z}^{(k,i)} := \text{POD}(\{\tilde{q}_j : \hat{I}_j^{c,(k)} = i\}, n, \|\cdot\|)$ for $i=1, \dots, k$.

4: **if** $\|\hat{I}^{c,(k)} - \hat{I}^{c,(k-1)}\| = 0$ **then**

5: **break**

6: **end if**

7: **end for**

References

[1] MMG Platform - Upgrade your meshes. https://www.mmgtools.org/, 2024.

[2] M. Aharon, M. Elad, and A. Bruckstein. K-SVD: An algorithm for designing overcomplete dictionaries for sparse representation. *IEEE Transactions on signal processing*, 54(11):4311–4322, 2006.

[3] D. Amsallem and B. Haasdonk. PEBL-ROM: Projection-error based local reduced-order models. *Advanced Modeling and Simulation in Engineering Sciences*, 3:1–25, 2016.

[4] D. Amsallem, M. J. Zahr, and C. Farhat. Nonlinear model order reduction based on local reduced-order bases. *International Journal for Numerical Methods in Engineering*, 92(10):891–916, 2012.

[5] D. Arthur and S. Vassilvitskii. K-means++: The advantages of careful seeding. volume 8, pages 1027–1035, 01 2007.

[6] O. Balabanov and A. Nouy. Randomized linear algebra for model reduction—part II: minimal residual methods and dictionary-based approximation. *Advances in Computational Mathematics*, 47:1–54, 2021.

[7] J. Barnett and C. Farhat. Quadratic approximation manifold for mitigating the Kolmogorov barrier in nonlinear projection-based model order reduction. *Journal of Computational Physics*, 464:111348, 2022.
[8] J. Barnett, C. Farhat, and Y. Maday. Neural-network-augmented projection-based model order reduction for mitigating the Kolmogorov bar-rier to reducibility. *Journal of Computational Physics*, 492:112420, 2023.

[9] N. Barral, T. Taddei, and I. Tifouti. Registration-based model reduction of parameterized PDEs with spatio-parameter adaptivity. *Journal of Computational Physics*, 499:112727, 2024.

[10] A. Bemporad and G. Cimini. Variable elimination in model predictive control based on K-SVD and QR factorization. *IEEE Transactions on Automatic Control*, 68(2):782-797, 2021.

[11] Y. Chen, T. T. Georgiou, and A. Tannenbaum. Optimal transport for Gaussian mixture models. *IEEE Access*, 7:6269-6278, 2018.

[12] D. S. Ching, P. J. Blonigan, F. Rizzi, and J. A. Fike. Reduced order modeling of hypersonic aerodynamics with grid tailoring. In *AIAA SCITECH 2022 Forum*, page 1247, 2022.

[13] T. Daniel, F. Casenave, N. Akkari, A. Ketata, and D. Ryckelynck. Physics-informed cluster analysis and a priori efficiency criterion for the construction of local reduced-order bases. *Journal of Computational Physics*, 458:111120, 2022.

[14] J. Delon and A. Desolneux. A Wasserstein-type distance in the space of Gaussian mixture models. *SIAM Journal on Imaging Sciences*, 13(2):936-970, 2020.

[15] F. Ducros, V. Ferrand, F. Nicoud, C. Weber, D. Darracq, C. Gacherieu, and T. Poinsot. Large-eddy simulation of the shock/turbulence interaction. *Journal of Computational Physics*, 152(2):517-549, 1999.

[16] R. Dupuis, J.-C. Jouhaud, and P. Sagaut. Surrogate modeling of aerodynamic simulations for multiple operating conditions using machine learning. *Aiaa Journal*, 56(9):3622-3635, 2018.

[17] J. L. Eftang, A. T. Patera, and E. M. Rønquist. An" hp" certified reduced basis method for parametrized elliptic partial differential equations. *SIAM Journal on Scientific Computing*, 32(6):3170-3200, 2010.

[18] D. Guignard and O. Mula. Tree-based nonlinear reduced modeling. *arXiv preprint arXiv:2404.12262*, 2024.

[19] T. Hastie, R. Tibshirani, and J. Friedman. *The Elements of Statistical Learning*. Springer Series in Statistics. Springer New York Inc., New York, NY, USA, 2009.

[20] A. Iollo and D. Lombardi. Advection modes by optimal mass transfer. *Physical Review E*, 89(2):022923, 2014.
[21] A. Iollo and T. Taddei. Mapping of coherent structures in parameterized flows by learning optimal transportation with Gaussian models. *Journal of Computational Physics*, 471:111671, 2022.

[22] L. Kaufman and P. J. Rousseeuw. *Finding groups in data: an introduction to cluster analysis*. John Wiley & Sons, 2009.

[23] Y. Kim, Y. Choi, D. Widemann, and T. Zohdi. A fast and accurate physics-informed neural network reduced order model with shallow masked autoencoder. *Journal of Computational Physics*, 451:110841, 2022.

[24] K. Lee and K. T. Carlberg. Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders. *Journal of Computational Physics*, 404:108973, 2020.

[25] M. A. Mirhoseini and M. J. Zahr. Model reduction of convection-dominated partial differential equations via optimization-based implicit feature tracking. *Journal of Computational Physics*, 473:111739, 2023.

[26] M. Ohlberger and S. Rave. Nonlinear reduced basis approximation of parameterized evolution equations via the method of freezing. *Comptes Rendus Mathematique*, 351(23-24):901-906, 2013.

[27] M. Ohlberger and S. Rave. Reduced basis methods: Success, limitations and future challenges. In *Proceedings of the Conference Algoritmy*, 2016. 1-12.

[28] A. Pinkus. *N-widths in Approximation Theory*, volume 7. Springer Science & Business Media, 2012.

[29] A. H. Razavi and M. Yano. Nonlinear model reduction for rapid and reliable CFD: application to transonic RANS flows. 2024.

[30] J. Reiss. Optimization-based modal decomposition for systems with multiple transports. *SIAM Journal on Scientific Computing*, 43(3):A2079-A2101, 2021.

[31] J. Reiss, P. Schulze, J. Sesterhenn, and V. Mehrmann. The shifted proper orthogonal decomposition: A mode decomposition for multiple transport phenomena. *SIAM Journal on Scientific Computing*, 40(3):A1322-A1344, 2018.

[32] T. Taddei. A registration method for model order reduction: data compression and geometry reduction. *SIAM Journal on Scientific Computing*, 42(2):A997-A1027, 2020.

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