# Error-Based mesh selection for efficient numerical simulations with variable parameters

Hugo Dornier, Olivier P Le Maître, Pietro Marco Congedo, Itham Salah El

Din, Julien Marty, Sébastien Bourasseau

► To cite this version:

Hugo Dornier, Olivier P Le Maître, Pietro Marco Congedo, Itham Salah El Din, Julien Marty, et al.. Error-Based mesh selection for efficient numerical simulations with variable parameters. 2025. [hal-05399931]

HAL Id: hal-05399931

https://hal.science/hal-05399931v1

Preprint submitted on 8 Dec 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
# Error-Based mesh selection for efficient numerical simulations with variable parameters

Hugo Dornier<sup>a,b,*</sup>, Olivier P. Le Maître<sup>c</sup>, Pietro M. Congedo<sup>b</sup>, Itham Salah el Din<sup>a</sup>, Julien Marty<sup>a</sup>, Sébastien Bourasseau<sup>d</sup>

<sup>a</sup>DAAA, ONERA, Institut Polytechnique de Paris, 92190 Meudon, France

<sup>b</sup>Inria, Centre de Mathématiques Appliquées, Ecole Polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

<sup>c</sup>CNRS, Centre de Mathématiques Appliquées, Ecole Polytechnique, IPP, Route de Saclay, 91128 Palaiseau Cedex, France

<sup>d</sup>DAAA, ONERA, Institut Polytechnique de Paris, 92322 Châtillon, France

## Abstract

Advanced numerical simulations often depend on mesh refinement techniques to manage discretization errors in complex models and reduce computational costs. This work concentrates on Adaptive Mesh Refinement (AMR) for steady-state solutions, which uses error estimators to iteratively refine the mesh locally and gradually tailor it to the solution. AMR requires evaluating the solution across a series of meshes. When solving the model for multiple operating conditions, such as in uncertainty quantification studies, full systematic adaptation can cause significant computational overhead. To mitigate this, the Error-based Mesh Selection (EMS) method is introduced to decrease the cost of adaptation. For each operating condition, EMS seeks to choose, from a library of pre-adapted meshes, the one that minimizes the discretization error. A key feature of this approach is the use of Gaussian Process models to predict the solution errors for each mesh in the library. These error models are built solely from the library's meshes and their solutions, using restriction errors as proxies for discretization errors, thereby avoiding additional model evaluations. The EMS method is tested on an analytical shock problem and a supersonic scramjet configuration, showing near-optimal mesh selection. The influence of library size on the resulting error level is also examined.

## 1. Introduction

Numerical simulations based on nonlinear physical models aim to provide accurate predictions by solving discretized governing equations. For steady problems, the quality of the numerical solution primarily depends on the spatial discretization of the computational domain. Resolving small-scale solution features typically requires locally refined meshes, while uniform refinement is computationally prohibitive. Mesh adaptation is therefore essential, yet the lack of a priori knowledge of the solution structure prevents prescribing an appropriate mesh beforehand.

Adaptive Mesh Refinement (AMR [1, 2, 3]) addresses this issue by iteratively adjusting the mesh resolution according to solution scales estimated on successive meshes. Combined with robust error estimators and refinement criteria, AMR enables the construction of meshes that meet accuracy requirements at controlled computational cost, and has been successfully applied to a wide range of complex problems [4, 5]. However, AMR incurs additional cost due to repeated solves, error evaluations, and mesh generation steps. This overhead becomes substantial in parametric, optimization, or sensitivity studies where many conditions must be evaluated [6, 7, 8]. To alleviate this burden, one would ideally employ a single mesh across all conditions. A common approach is to adapt the mesh to a nominal condition, but when the solution is highly parameter-dependent, this may lead to large errors compared with systematic adaptation [9, 10].

<footnote>*Corresponding author
Email address: hugo.dornier@gmail.com (Hugo Dornier)</footnote>
Adjoint-based enrichments [11] partially extend the validity of nominal meshes but remain fundamentally local and cannot guarantee error control when solution structures vary markedly with the parameters.

Alternative strategies aim at constructing a single mesh representative of several conditions. In [12], local mesh metrics obtained from independently adapted meshes were combined to generate a unified mesh. This approach is more robust than nominal adaptation, though its performance depends on the choice of construction conditions. More recently, Mean Mesh Adaptation (MMA [13]) was introduced to control the average error over a continuous parameter set via Monte Carlo estimates. MMA performs well when refinement remains confined to limited regions of the domain; however, when solution features shift widely across the domain, the method yields near-uniform meshes and becomes less efficient than systematic AMR. Overall, methods relying on a single mesh are inherently restricted to problems with relatively stable solution structures.

The present work adopts a different perspective by considering multiple meshes, each adapted to a specific condition. These adapted meshes, together with their solutions and parameter values, form a library assumed to be given. We propose an Error-based Mesh Selection (EMS) method that, for any new condition, selects a priori the mesh from this library expected to yield the smallest error.

EMS is inspired by quantization theory [14, 15], which maps high-dimensional input data onto a finite codebook. Specifically, EMS parallels vector quantization techniques that partition a domain by minimizing a metric, as in clustering methods (e.g., k-means). Here, the mapping assigns each condition in a continuous parameter space to one mesh in the library while minimizing the associated discretization error. Assuming the availability of projection and interpolation operators between the approximation spaces, as well as a computable a posteriori error estimator, EMS constructs error models for each mesh directly from the library. A key advantage is that no additional flow solutions are required during the selection phase.

The remainder of the paper is organized as follows. Section 2 reviews the discretization framework, the class of problems considered, and the mesh adaptation strategies relevant to this work. Section 3 introduces the EMS formulation and the proposed approximation strategy. Section 4 presents an analytical CFD test case involving a supersonic oblique shock to illustrate the method's performance over a wide parameter range. Section 5 extends the study to a realistic supersonic scramjet configuration [16], demonstrating the robustness and practical relevance of the approach.

## 2. Background materials

The EMS method proposed below is largely independent of the specific numerical schemes and adaptation strategies employed to obtain the discrete solutions and generate the adapted meshes. In this section, we state the general requirements for these methods to be compatible with EMS and introduce the notations used throughout the article.

### 2.1. Mesh and discrete solution

We consider a fixed domain $\Omega \subset \mathbb{R}^d$, $d \ge 1$, with boundary $\partial\Omega$. Let $u(x) \in \mathcal{V}(\Omega)$ denote the unique solution of the problem, which satisfies the equations

$$ \mathcal{P}(u) = 0, \qquad (1) $$

Here, the operator $\mathcal{P}$ encompasses both the governing equations and the boundary conditions. The solution space $\mathcal{V}(\Omega)$ depends on the particular problem and is assumed to be equipped with a norm $\| \cdot \|_{\mathcal{V}}$. In what follows, the dependence of $\mathcal{V}$ on $\Omega$ is omitted for convenience.

To solve (1), we employ a discretization of $\Omega$ represented by a mesh $\mathcal{H}$. Without loss of generality, $\mathcal{H}$ is defined as a collection of non-overlapping computational elements $c_i$ that partition $\Omega$:

$$ \Omega = \bigcup_{i \in \mathcal{H}} c_i, \quad \forall (c_i, c_j) \in \mathcal{H}^2 c_i \cap c_j = \emptyset \text{ for } i \neq j. $$

This definition encompasses, for instance, the finite-volume, finite-element, and discontinuous Galerkin methods. As is standard in the literature, $h(\mathcal{H})$ denotes the size of the largest computational element,
$h(\mathcal{H}) = \max_{c_i \in \mathcal{H}} |c_i|$, while $|\mathcal{H}|$ denotes the number of computational elements in the mesh. We introduce the approximation space $V_{\mathcal{H}} \subset V(\Omega)$ constructed on the mesh, and denote by $u_{\mathcal{H}} \in V_{\mathcal{H}}$ the discrete solution. The dimension of the approximation space, $\dim(V_{\mathcal{H}})$, is typically larger than $|\mathcal{H}|$, even for scalar problems.

We assume that the discrete problem is well-posed with a unique solution, and that $\|u_{\mathcal{H}} - u\|_{V} \to 0$ as $h(\mathcal{H}) \to 0$.

In practice, it is neither possible nor desirable to let $h(\mathcal{H})$ tend to zero to achieve a sufficiently accurate solution. A more effective strategy is to control the size of individual computational elements to attain a prescribed accuracy with minimal computational effort. For the discussion that follows, we assume that the objective is to control the approximation error measured in the $L_2(\Omega)$ norm:

$$ \varepsilon_{\mathcal{H}}^2 = \|u - u_{\mathcal{H}}\|_{L_2(\Omega)}^2 = \sum_{c_i \in \mathcal{H}} \|u - u_{\mathcal{H}}\|_{L_2(c_i)}^2. \qquad (2) $$

However, error control based on other error definitions may also be considered. For instance, in goal-oriented adaptation, one may wish to control the error,

$$ \varepsilon_{\mathcal{H}} = |\Phi(u) - \Phi(u_{\mathcal{H}})|, $$

where $\Phi: V \to \mathbb{R}$ is a quantity of interest.

## 2.2. Adaptive Mesh Refinement (AMR)

Mesh adaptation aims to determine a mesh $\mathcal{H}$ that (approximately) minimizes the mesh size $|\mathcal{H}|$ while satisfying the error condition $\varepsilon_{\mathcal{H}} \le \varepsilon_{ad}$, where $\varepsilon_{ad} > 0$ is a prescribed target error. Here, we assume that the computational cost of obtaining $u_{\mathcal{H}}$ increases monotonically with $\dim(V_{\mathcal{H}})$, which is proportional to $|\mathcal{H}|$.

In AMR methods, a sequence of refined meshes $\mathcal{H}^{(0)}, \mathcal{H}^{(1)}, \dots$ is constructed until the error criterion is satisfied. In the simplest algorithms, the current mesh $\mathcal{H}^{(k)}$ is refined to $\mathcal{H}^{(k+1)}$ by subdividing selected computational elements, so that $V_{\mathcal{H}^{(k)}} \subset V_{\mathcal{H}^{(k+1)}}$. A criterion is required to prioritize the elements that need refinement and to avoid refining elements $c_i$ that already contribute little to the error. A common choice is to select for refinement a subset comprising a fraction $\alpha_{ad}$ of the computational elements in $\mathcal{H}^{(k)}$, with $0 < \alpha_{ad} < 1$, based on the local error contributions.

For instance, in the case of the $L_2$ error norm in (2) and a scalar solution $u$, the decomposition reads

$$ \varepsilon_{\mathcal{H}^{(k)}}^2 = \sum_{c_i \in \mathcal{H}^{(k)}} \varepsilon_{\mathcal{H}^{(k)}}^{\text{loc}}(c_i)^2, \quad \varepsilon_{\mathcal{H}^{(k)}}^{\text{loc}}(c_i)^2 = \|u - u_{\mathcal{H}^{(k)}}\|_{L_2(c_i)}^2 = \int_{c_i} |u - u_{\mathcal{H}^{(k)}}|^2 dx. $$

In practice, the exact solution $u$ is unknown, and estimators $\hat{\varepsilon}_{\mathcal{H}}$ of $\varepsilon_{\mathcal{H}}$ are introduced, for instance using a posteriori error estimation [17, 18, 19, 2]. In the case of goal-oriented adaptation, an additional adjoint problem may need to be solved to estimate the contribution of each computational element to the error.

Once the fraction $\alpha_{ad}$ of computational elements carrying the highest estimated error $\hat{\varepsilon}_{\mathcal{H}^{(k)}}$ is determined, the next mesh $\mathcal{H}^{(k+1)}$ is generated through a refinement procedure. Various types of refinement can be employed. One approach is to subdivide the selected elements $c_i$, possibly including some neighboring elements to maintain a smooth variation of cell sizes $|c_i|$ and ensure mesh quality. More generally, isotropic or anisotropic metrics can be defined to constrain the local sizes of the computational elements in the refined mesh $\mathcal{H}^{(k+1)}$, in which case the nested character of the meshes and approximation spaces may no longer hold.

During the iterative process, the parameter $\alpha_{ad}$ plays a crucial role. For small values of $\alpha_{ad}$, only a few elements are refined at each iteration, so that $\mathcal{H}^{(k)}$ and $\mathcal{H}^{(k+1)}$ differ only slightly. Many iterations may then be required to obtain a mesh satisfying the target error $\varepsilon_{ad}$, but the resulting mesh will be sharply adapted. In contrast, for larger $\alpha_{ad}$, more elements are refined per iteration, which may yield a mesh that is less sharply adapted, though it can be reached in fewer iterations. In the present work, we assume that $\alpha_{ad}$ is chosen appropriately to ensure well-adapted meshes and stays fixed throughout adaptation.

The structure of the AMR algorithm is summarized in Algorithm 1.
**Algorithm 1** AMR algorithm.

**Require:** $\mathcal{H}^{(0)}$, $\varepsilon_{ad}$, $\alpha_{ad}$

1: Set $k = 0$

2: Compute $u_{\mathcal{H}^{(k)}}$ and local error estimates $\hat{\varepsilon}_{\mathcal{H}^{(k)}}^{\text{loc}}(c_i)$

3: **while** $\sum_{c_i \in \mathcal{H}^{(k)}} \hat{\varepsilon}_{\mathcal{H}^{(k)}}^{\text{loc}}(c_i)^2 > \varepsilon_{\text{ad}}^2$ **do**

4: Flag the $\alpha_{\text{ad}}|\mathcal{H}^{(k)}|$ elements with highest error estimates

5: Refine the flagged elements of $\mathcal{H}^{(k)}$ to generate $\mathcal{H}^{(k+1)}$

6: $k \leftarrow k + 1$

7: Compute $u_{\mathcal{H}^{(k)}}$ and local error estimates $\hat{\varepsilon}_{\mathcal{H}^{(k)}}^{\text{loc}}(c_i)$

**return** $\mathcal{H}^{(k)}, u_{\mathcal{H}^{(k)}}$

## 2.3. Present methods

As outlined in the introduction to this section, EMS is independent of the specific numerical methods employed for problem resolution and mesh refinement. Nevertheless, specific choices were made for its implementation and testing in this work. To investigate the performance of the EMS method, we consider CFD problems, and more specifically, supersonic non-viscous configurations. These configurations are particularly challenging for studies involving variable conditions, as they feature local phenomena such as shocks and expansion fans that are highly sensitive to boundary conditions. Consequently, this study employs a standard Finite Volume Method (FVM) with a first-order Roe scheme [20], which ensures stable handling of solution discontinuities such as shocks while maintaining a satisfactory level of accuracy.

The local error estimator guiding the adaptation process is based on an *a posteriori* estimation of the interpolation error on $u_{\mathcal{H}}$. The specific estimator used here is a piecewise constant interpolation error estimator derived from our previous work [13], tailored for piecewise constant schemes such as the Roe scheme employed in this study. The error in a given cell $c_i$ is defined in Equations 3 and 4 of [13]. More generally, for $n$-th order schemes, the interpolation error estimate relies on the $(n+1)$-th derivative of the solution.

For adaptation, we adopt a standard isotropic subdivision strategy, in which cells marked for refinement are subdivided into $2^d$ sub-cells. A smoothing procedure is then applied to the resulting mesh to ensure a gradual evolution of cell sizes, so that the difference in refinement level between two neighboring cells does not exceed 1. This approach produces a sequence of nested meshes $\mathcal{H}^{(k)}$ starting from an initial coarse mesh $\mathcal{H}^{(0)}$, which can be conveniently represented using a tree structure. This structure can be efficiently exploited for compact adaptation procedures, storage, and various operations such as initial solution transfer, similarity measurements, and multi-resolution computations. Finally, we note that mesh coarsening is not considered in this study.

## 3. Error-based mesh selection

In this section, we extend the problem of predicting $u$ to multiple conditions. We introduce a vector of parameters $q \in \mathcal{Q} \subseteq \mathbb{R}^{d_q}$ that defines these conditions. For example, $q$ may parametrize boundary conditions, model constants, or external forcings. Let $\mu_{\mathcal{Q}}$ denote the probability measure of $q$, such that

$$ \mathbb{E}_q [f(q)] = \int_{\mathcal{Q}} f(q) d\mu_{\mathcal{Q}}(q), $$

for any quantity of interest $f(q)$.

We restrict ourselves to a fixed computational domain $\Omega$ (independent of $q$) and denote by $u(q) \in \mathcal{V}$ the solution satisfying the parameter-dependent problem

$$ \mathcal{P}(u, q) = 0. \qquad (3) $$

We assume that (3) is almost surely well-posed, and denote by $u_{\mathcal{H}}(q) \in \mathcal{V}$ and $\varepsilon_{\mathcal{H}}(q)$ the discrete solution and associated error corresponding to its resolution on a mesh $\mathcal{H}$.
### 3.1. Library of meshes

In the following, $\mathcal{H}_q$ denotes the mesh adapted to $u(q)$, i.e., to the resolution of (3) for a given $q$. The mesh $\mathcal{H}_q$ is constructed such that the estimated error

$$ \|u(q) - u_{\mathcal{H}(q)}(q)\| \approx \hat{\varepsilon}_{\mathcal{H}(q)}(q) \leq \varepsilon_{\text{ad}}. \quad (4) $$

Instead of computing $\mathcal{H}_q$ for each queried parameter $q \in \mathcal{Q}$, or using a single mesh $\bar{\mathcal{H}}$ for all $q$, we propose to address the problem by selecting a mesh from a finite library of $1 < N_L < \infty$ meshes.

Let $\{q_i, i = 1, \dots, N_L\}$ be a set of distinct parameter vector values in $\mathcal{Q}$. For each $i = 1, \dots, N_L$, we define $\mathcal{H}_i \doteq \mathcal{H}(q_i)$ and $u_i \doteq u_{\mathcal{H}_i}(q_i)$. The task considered in this section is then to select the appropriate mesh $\mathcal{H}_i$ from the library to approximate the solution $u(q)$ for a queried parameter $q \in \mathcal{Q}$.

Assuming that a mesh $\mathcal{H}_i$ is suitable in the neighborhood of $q_i$, a simple approach for mesh selection relies on the parametric distance. Specifically, we choose the mesh $\mathcal{H}_{i_{\text{dist}}}$ with index $i_{\text{dist}}$ given by

$$ i_{\text{dist}} \in \arg \min_{i \in \{1, \dots, N_L\}} d_W(q_i, q), \quad (5) $$

where $d_W$ is the weighted Euclidean distance on $\mathbb{R}^{d_q}$,

$$ d_W(q, q')^2 = (q - q')^T W (q - q'), \quad (6) $$

where $W$ is a symmetric positive definite matrix. Since the components of the vector $q$ may correspond to different physical quantities, the matrix $W$ is essential for properly scaling these quantities. Transformations that map the components' ranges to the unit interval do not reduce $W$ to the identity matrix, as the distance would still depend on the parametric ranges. In practice, $W$ should be fitted to the specific problem under consideration, a task that can require considerable effort and may yield limited efficiency. This observation motivates the development of a more robust approach in the remainder of this section.

### 3.2. Error-based mesh selection

We formulate the problem in terms of a policy $\mathfrak{p}: \mathcal{Q} \to \{1, \dots, N_L\}$ that aims to minimize the resulting approximation error

$$ \varepsilon_{\mathcal{H}_{\mathfrak{p}(q)}}(q) = \|u(q) - u_{\mathcal{H}_{\mathfrak{p}(q)}}(q)\|. \quad (7) $$

The solution $u_{\mathcal{H}_{\mathfrak{p}(q)}}(q)$ represents the approximate solution of the problem with parameter $q$, computed on the mesh from the library selected by the policy at $q$.

To determine the optimal policy, we must specify the class of policies under consideration and clarify the sense in which optimality is defined. Since $q$ is a random variable and the error is non-negative, any moment of its distribution could serve as a valid criterion. A natural choice is to minimize the mean error for the policy:

$$ \mathcal{E}(\mathfrak{p}) = \mathbb{E}_q [\varepsilon_{\mathcal{H}_{\mathfrak{p}(q)}}(q)]. \quad (8) $$

However, evaluating the performance of a policy typically requires computing the error for many values of $q$, which makes this approach impractical.

To derive an alternative approach, we observe that for any policy $\mathfrak{p}$, a lower bound on the mean error in (8) is

$$ \mathcal{E}(\mathfrak{p}) \geq \mathbb{E}_q \left[ \min_{i=1, \dots, N_L} \varepsilon_{\mathcal{H}_i}(q) \right], \quad \varepsilon_{\mathcal{H}_i}(q) = \|u(q) - u_{\mathcal{H}_i}(q)\|. $$

This expression highlights the role of the error functions $\varepsilon_{\mathcal{H}_i}(q)$ in determining the performance of a policy. Rather than directly modeling $\mathfrak{p}$, we propose to model selected error functions to construct an approximate policy. The motivation is that these error functions can be inferred from the library without any additional problem resolution, as explained below and illustrated in the results section.

Computing the error functions for each library mesh over all queried $q$ is not feasible, since $u(q)$ is unknown. The choice is made instead to rely on the error estimates $\hat{\varepsilon}_{\mathcal{H}_i}(q)$, which would require solving
$u_{\mathcal{H}_i}(q)$ on the library meshes. Performing these $N_L$ resolutions may not be practical compared to performing a full specific adaptation, unless $N_L$ is small. To overcome this limitation, we introduce a proxy error model for $\varepsilon_{\mathcal{H}_i}(q)$, which can be determined directly from the library.

### 3.3. Error models

A first approach to constructing the error models from the library is to compute, for all $1 \le i, j \le N_L$, the estimated errors $e_{i,j} \doteq \hat{\varepsilon}_{\mathcal{H}_i}(q_j)$, and use them to build a model $e_i(\mathbf{p})$. This approach requires $N_L(N_L - 1)$ additional resolutions on known meshes. To reduce this computational effort, we further approximate $\hat{\varepsilon}_{\mathcal{H}_i}(q_j)$ using the so-called restriction error [21] formulation.

We define $\Pi_i$ as the projection operator from $\mathcal{V}$ onto $\mathcal{V}_{\mathcal{H}_i}$, given by

$$ \Pi_i u = \arg \min_{v \in \mathcal{V}_{\mathcal{H}_i}} \|u - v\|. \qquad (9) $$

This defines an orthogonal projection since $\mathcal{V}_{\mathcal{H}_i} \subset \mathcal{V}$. This projection can be computed efficiently in the case of MRA (multi-resolution analysis) methods, particularly when the meshes in the library are structured and represented by (possibly incomplete) trees. In fact, most MRA libraries provide utilities to project functions from one mesh to another. Writing the error $\varepsilon_{\mathcal{H}_i}(q)$ at $q = q_j$, we obtain

$$ \begin{aligned} \varepsilon_{\mathcal{H}_i}(q_j)^2 &= \|u(q_j) - u_{\mathcal{H}_i}(q_j)\|^2 \le \|u(q_j) - u_j\|^2 + \|u_j - u_{\mathcal{H}_i}(q_j)\|^2 \\ &\le \varepsilon_{\text{ad}} + \|u_j - \Pi_i u_j\|^2 + \|\Pi_i u_j - u_{\mathcal{H}_i}(q_j)\|^2. \end{aligned} \qquad (10) $$

The second term, $\|u_j - \Pi_i u_j\|$, is known as the restriction error and is denoted by $e_i$ in the following. It is defined such that the restriction error vanishes when $\mathcal{V}_{\mathcal{H}_i}$ is locally richer (more refined) than $\mathcal{V}_{\mathcal{H}_j}$. Consequently, it quantifies the loss of quality in $u(q_j)$ when projected onto $\mathcal{V}_{\mathcal{H}_i}$. This quantity can be computed with limited computational effort, since the adapted solutions $u_j$ are readily available.

The final term, $\|\Pi_i u_j - u_{\mathcal{H}_i}(q_j)\|$, is generally expected to be dominated by the restriction error, as it quantifies the difference between a reasonably accurate approximation of the projection of $u(q_j)$ onto $\mathcal{V}_{\mathcal{H}_i}$ and the actual solution computed on $\mathcal{H}_i$.

The next task is to construct surrogate models $\hat{e}_i(q)$ for the restriction errors. These restriction errors can be efficiently computed from the library for $q = q_j$, $j = 1, \dots, N_L$, as explained above, providing the data necessary to train the models. Here, we employ Gaussian process (GP) regression [22] to approximate the restriction error functions. This choice is motivated by the robustness of GP models, which can accommodate limited training data. Additionally, GP models provide confidence intervals that can be leveraged to develop infilling strategies for enhancing the library in future work.

The models take the following generic form:

$$ e_i(q) \approx \hat{e}_i(q) = \sum_{k=0}^{n} \beta_k \hat{f}_k(q) + Z(q). \qquad (11) $$

In (11), the $\hat{f}_k$ are prescribed trend functions, and $Z$ is a Gaussian process. The GP, $Z$, has a prior distribution, which is updated via Bayes' formula using the $N_L$ observations $e_i(q_j)$. The prior distribution consists of a mean function and a covariance function $C(q, q')$, generally involving hyperparameters that must be set alongside the trend coefficients $\beta_k$. In the present work, we rely on the library smt [23], which provides several choices of trend functions and covariance structures, as well as different algorithms for regression and hyperparameter optimization. Based on numerical tests, we train the GP models on standardized data (centered and scaled), using linear trend functions, zero-mean priors, and simple squared exponential covariance functions.

$$ C(q, q') = \sigma^2 \exp \left[ - \frac{\|q - q'\|^2}{l_c^2} \right], \qquad (12) $$

with variance ($\sigma^2$) and correlation length ($l_c$) as hyperparameters.
One challenge in approximating the restriction error $e_i$ from a limited dataset $e_i(q_j)$ is that the approximation must be sufficiently accurate when the error is small, to discriminate between meshes that induce low errors. For a large library, the restriction error $e_i$ may be small only in a limited neighborhood of $q_i$ and large elsewhere, so that most data points $e_i(q_j)$ take large values. Training a GP on a dataset dominated by high-error points can compromise the quality of the prediction. To mitigate this effect, we found it more effective to train the GP model $\hat{e}_i(q)$ using only the restriction errors corresponding to parameters $q_j$ within a fixed neighborhood $N_i$ of $q_i$. Defining this neighborhood geometrically based on the distance $\|q_i - q_j\|$ encounters the same difficulties discussed in 3.1, and we adopted the following definition,

$$N_i = \{j = 1, \dots, N_L | e_i(q_j) \le \delta_e\}, \quad (13)$$

where $\delta_e$ is a user-prescribed error threshold. With this definition, the GP training disregards data points with excessively large errors, thereby improving the predictive capability in regions where the error is small.

### 3.4. Approximate EMS method

Given a library of meshes (and their associated solutions), the approximate EMS (A-EMS) method proceeds in two stages. In the offline stage, the restriction error models $\hat{e}_i(q)$ are constructed. In the in-line stage, the index of the mesh to use for a queried parameter value $q$ is selected from

$$\mathfrak{p}(q) = \arg \min_{i=1, \dots, N_L} \hat{e}_i(q).$$

The selected mesh will therefore be the one with the lowest restriction error predicted by the GP.

To assess the performance of the A-EMS method, we focus on the resulting average error

$$\mathcal{E}(\mathfrak{p}) = \mathbb{E}_q [\varepsilon_{\mathcal{H}_{\mathfrak{p}(q)}}(q)]. \quad (14)$$

The A-EMS error can be compared to the average minimal error achievable with the considered library. Let $i^*(q)$ denote the true optimal mesh selection,

$$i^*(q) = \arg \min_{i=1, \dots, N_L} \|u(q) - u_{\mathcal{H}_i(q)}\|, \quad (15)$$

and denote corresponding $\varepsilon^*(q)$ the minimal error,

$$\varepsilon^*(q) = \|u(q) - u_{\mathcal{H}_{i^*(q)}}(q)\|.$$

With these notations, the local sub-optimality of the A-EMS method is defined by

$$\rho_{\text{A-EMS}}(q) = \frac{\varepsilon_{\mathcal{H}_{\mathfrak{p}(q)}}(q) - \varepsilon^*(q)}{\varepsilon_{\text{ad}}}, \quad (16)$$

and the global sub-optimality is

$$\bar{\rho}_{\text{A-EMS}} = \frac{\|\varepsilon_{\mathcal{H}_{\mathfrak{p}(q)}}(q) - \varepsilon^*(q)\|^2}{\varepsilon_{\text{ad}}}. \quad (17)$$

To complete the assessment of the A-EMS method, we provide a brief analysis of its computational cost and compare it with the cost of performing a systematic mesh adaptation for all queried parameters $q$, referred to as Specific Mesh Adaptation (SMA). The computational time required to generate the library in this analysis is not taken into account, as it is subsequently amortized over a large number of queried values $q$, which is assumed to be much greater than $N_L$.

Under the working assumption of an MRA method, where the mesh $\mathcal{H}^{(k+1)}$ is obtained by refining a fixed fraction $\alpha_{ad}$ of the cells in $\mathcal{H}^{(k)}$ into $2^d$ sub-cells, we obtain

$$|\mathcal{H}^{(k)}| = |\mathcal{H}^{(0)}| G_{d, \alpha_{ad}}^k, \quad G_{d, \alpha_{ad}} = (1 + (2^d - 1)\alpha_{ad})$$
where $|\mathcal{H}^{(0)}|$ is the number of cells in the initial mesh, $G_{d,\alpha_{ad}} > 1$ is the mesh size increase factor per adaptive iteration, $d$ is the number of spatial dimensions, and $\alpha_{ad}$ is the fraction of cells refined at each iteration.

To estimate the computational cost of the methods, we assume that the cost of error estimation and mesh refinement is negligible compared to solving the problem on a given mesh, which is often the case with modern (possibly parallel) mesh adaptation frameworks. Therefore, for systematic mesh adaptation, the total cost to obtain a solution on a sufficiently refined mesh corresponds to the cost of solving the problems on the sequence of refined meshes. We further assume that the solution time for a given mesh can be expressed as

$$T(\mathcal{H}) = C|\mathcal{H}|^{\gamma},$$

with $\gamma \ge 1$ characterizing the solver's scalability. The total computational time to obtain the solution on the adapted mesh $\mathcal{H}^{(k_{ad})}$ is then

$$T_{k_{ad}} = C|\mathcal{H}^{(0)}|^{\gamma} \sum_{k=0}^{k_{ad}} (G_{d,\alpha_{ad}}^{k})^{\gamma},$$

where $k_{ad}$ is the number of adaptation iterations required to reach the prescribed accuracy. In practice, $k_{ad}$ depends on the parameter $q$, the target precision $\varepsilon_{ad}$, and the refinement fraction $\alpha_{ad}$. We are interested in the average computational time to solution, $\bar{T}_{SMA} = \mathbb{E}_q [T_{k_{ad}}]$. For small $k_{ad}$ variations with $q$, the SMA average time to solution can be approximated by

$$\bar{T}_{SMA} \approx C|\mathcal{H}^{(0)}|^{\gamma} \sum_{k=0}^{\bar{k}_{ad}} (G_{d,\alpha_{ad}}^{k})^{\gamma}, \quad (18)$$

where $\bar{k}_{ad}$ is the (rounded) average number of adaptive iterations.

For the A-EMS method, only a single solution needs to be computed for each queried $q$, with a cost

$$T_{\text{A-EMS}}(q) = C|\mathcal{H}_{p(q)}|^{\gamma}.$$

When the library consists of $N_L$ samples drawn uniformly from $d\mu_Q$, we have $\mathbb{E}_q [|\mathcal{H}_i|] \approx |\mathcal{H}^{(0)}|G_{d,\alpha_{ad}}^{\bar{k}_{ad}}$. It is also reasonable to assume that $\mathbb{E}_q [|\mathcal{H}_{p(q)}|] \approx |\mathcal{H}^{(0)}|G_{d,\alpha_{ad}}^{\bar{k}_{ad}}$ if the policy performs correctly. Under these assumptions, the average time to solution for the A-EMS method is

$$\bar{T}_{\text{A-EMS}} \approx C|\mathcal{H}^{(0)}|^{\gamma} (G_{d,\alpha_{ad}}^{\bar{k}_{ad}})^{\gamma}. \quad (19)$$

The averaged gain in time per solution of the A-EMS method, compared to SMA, is

$$\bar{T}_{\text{SMA}} - \bar{T}_{\text{A-EMS}} \approx C|\mathcal{H}^{(0)}|^{\gamma} \sum_{k=0}^{\bar{k}_{ad}-1} (G_{d,\alpha_{ad}}^{k})^{\gamma} = C|\mathcal{H}^{(0)}|^{\gamma} \frac{1 - (G_{d,\alpha_{ad}}^{\gamma})^{\bar{k}_{ad}}}{1 - G_{d,\alpha_{ad}}^{\gamma}}.$$

This expression indicates that the differences will primarily depend on the refinement fraction $\alpha_{ad}$ and the average number of adaptation iterations $\bar{k}_{ad}$, two quantities that are not independent, since reducing $\alpha_{ad}$ generally requires more adaptation iterations. The reduction in computational time will also increase as the target precision $\varepsilon_{ad}$ decreases, because $\bar{k}_{ad}$ will increase. However, the average error of A-EMS does not necessarily decrease unless the library is sufficiently rich.

## 4. EMS investigation on oblique shock

This section examines the proposed selection method using an analytical problem that permits the exact projection of the solution onto a prescribed mesh and the evaluation of the resulting error at low computational cost.
### 4.1. Oblique shock problem presentation

We consider the two-dimensional flow of an inviscid, compressible fluid around a symmetric wedge of fixed angle $\theta_w$, see Figure 1a. For a supersonic inflow, an attached oblique shock forms symmetrically at the wedge tip and extends outward at a constant angle $\beta$, which depends on the inflow Mach number $M_0$, the wedge angle $\theta_w$, and the ratio of specific heats of the perfect gas, $\gamma$.

The oblique shock separates two constant regions: the inflow region and the post-shock region. Both are known explicitly as functions of the problem parameters ($\theta_w, M_0, \gamma$), see [24]. Figure 1b illustrates the dependence of the shock angle $\beta$ on $M_0$ and $\gamma$ for $\theta_w = 15^\circ$. The shock angle decreases with $M_0$ and increases with $\gamma$. The effect of $\gamma$ becomes more pronounced at higher Mach numbers: it is almost negligible around $M_0 \approx 2$, but dominant near $M_0 \approx 8$.

This variation in the dependence of $\beta$ on the flow parameters poses a challenge for mesh selection methods, since a simple parametric distance cannot capture such heterogeneous dependencies. For the analysis, we consider two ranges of uncertain conditions. The reference range is $(M_0, \gamma) \in R_2 = [2, 8] \times [1.4, 2]$, and the restricted range is $R_1 = [2, 4] \times [1.4, 1.6]$, both assumed with uniform probability distributions.

<figure><img src="image_2.png" /><figcaption>Figure 1: Presentation of the oblique shock problem.</figcaption></figure>

This simple flow configuration enables a cost-effective and accurate piecewise-constant approximation on any prescribed mesh, enabling precise evaluation of discretization errors. For mesh adaptation, we focus only on the Mach field and employ a standard AMR strategy that isotropically subdivides a fraction $\alpha_{ad}$ of the mesh cells with the largest local error estimates, continuing until the estimated $L_2$ error falls below $\varepsilon_{ad} = 0.005$. All computations presented use the interpolation error for error estimation. The refinement fraction is fixed at $\alpha_{ad} = 0.02$ to produce sharply adapted meshes.

### 4.2. Selection maps and error analysis

In this section, we consider a library of $N_L = 16$ meshes associated with parameter values $q_i$, generated by a Sobol sequence in $R_2$. Figure 2 reports the positions of the parameters $q_i$ using symbols, while subdomains of $R_2$ sharing the same minimizer $i^*$ in equation (15) are shown with a common color. We distinguish between used and unused meshes by the marker shape: meshes selected at least once are represented at their parameter location $q_i^*$ with hexagons, whereas meshes never selected are indicated with triangles. The three maps in Figure 2 are generated by evaluating the errors over an $80 \times 80$ uniform grid of parameter values.

The first plot, Figure 2a, shows the optimal selection map. A first observation is that some library meshes are never selected. This occurs because a more refined mesh, adapted to another condition, may yield a similar shock angle $\beta$. Since the jump across the shock increases with $M_0$, the optimal strategy favors the mesh adapted to the higher Mach number when two candidates produce comparable shock angles. Furthermore, because the $q_i$ are drawn randomly, the resulting selection domains exhibit highly disparate sizes, reflecting the nonlinear parametric behavior highlighted above.
Figure 2b presents the selection map produced by the A-EMS method. For comparison, the contours of the optimal map are also shown. One observes that, while the A-EMS map captures some features of the optimal map, such as the Mach-aligned anisotropy and the increasing influence of $\gamma$ as $M_0$ grows, the boundaries of the selection domains are generally misplaced. In addition, all but one mesh are selected, and the resulting domains are sometimes disconnected or small. The GP modeled errors are likely responsible for these complex selection domains, as they can lead to a noisy policy in areas where the restriction errors of different meshes are comparable, particularly at the boundaries of the optimal selection domains. Nevertheless, as discussed below, selecting a non-optimal mesh remains acceptable when the induced error is close to that of the optimal choice.

This relatively coherent behavior contrasts with the selection map of the distance-based mesh selection (DMS) method, shown in Figure 2c. Here, the map is obtained using a diagonal weight matrix derived from the normalization of the variation ranges of $M_0$ and $\gamma$. The resulting partition of $Q$ is simply the Voronoi diagram induced by $W$. In this case, all meshes are selected within their parametric neighborhood; thus, the domains’ boundaries differ substantially from those of the optimal map.

<figure><img src="image_3.png" /><figcaption>Figure 2: Selection maps for the oblique shock problem on $R_2$ with $N_L = 16$. The library meshes parameters $q_i$ are shown with symbols.</figcaption></figure>

While indicative of proximity to the optimal selection, the maps are not directly informative of the errors achieved by the A-EMS and DMS strategies. To complete the analysis, we examine the local errors for each method in Figure 3. The figure reports the local sub-optimality errors normalized by the target precision $\epsilon_{ad}$, defined in (16). Comparing Figures 3a and 3b, we observe that the sub-optimality of A-EMS remains below $\epsilon_{ad}$ across most of the domain, whereas for DMS it exceeds $2\epsilon_{ad}$ in several regions (note that different color scales are used for each). Moreover, the regions of high sub-optimality for A-EMS tend to concentrate along the boundaries of the selection domains, i.e., far from the library meshes. For both methods, the sub-optimality gradients are primarily aligned with the Mach direction, underscoring the greater mesh discrepancies along this parameter and the resulting reduction in accuracy of both modeling and selection in this direction.

Overall, the sub-optimality map of A-EMS shows that when a non-optimal mesh is selected, the induced error is typically not much larger than the optimal one, particularly in regions where the library density is sufficient. Notable exceptions occur in the bottom-left and top-right portions of $Q$, where large sub-optimality values are observed. In these areas, the policy based on the modeled restriction errors selects meshes adapted to significantly different shock angles. This behavior stems from the degradation of the GP approximations near domain boundaries due to extrapolation. It could be mitigated by sequentially enriching the mesh library.

Figure 3c further compares the error distributions obtained with the optimal selection, A-EMS, and DMS. Interestingly, the probability of achieving an error smaller than $\epsilon_{ad}$ is nonzero for all methods, particularly for
the optimal selection and A-EMS. This effect arises from overlapping shock angles, which lead to resolution on overadapted meshes. The A-EMS error distribution exhibits a much lighter tail than that of DMS, with an almost negligible probability of exceeding $4\epsilon_{ad}$. This result highlights the robustness of A-EMS, which is unlikely to produce excessively large errors.

<figure><img src="image_4.png" /><figcaption>Figure 3: Sub-optimality maps and error distributions. Oblique shock on $R_2$ with $N_L = 16$.</figcaption></figure>

To obtain a more quantitative comparison of the selection strategies, Table 1 reports the evolution of the averaged error for different library sizes and for the two parametric ranges. For convenience, these averaged errors are normalized by the target error $\epsilon_{ad}$. A first observation is that the averaged error decreases as $N_L$ increases and converges toward $\epsilon_{ad}$. As expected, the convergence of A-EMS is slightly slower than that of the optimal selection but significantly faster than that of DMS, with this trend being more pronounced for the larger range. Finally, the lower averaged error of the optimal selection for $R_2$ compared to $R_1$ at $N_L = 64$ is explained by the stronger effect of shock angle overlap in the wider range.

<table><caption>Table 1: Averaged error $\mathcal{E}(p)$, normalized by the target error $\epsilon_{ad}$ for the optimal selection (Opt.), A-EMS, and DMS methods. Using different library sizes $N_L$ and the two parametric ranges. Oblique shock problem.</caption><thead><tr><th colspan="10">×ε<sub>ad</sub></th></tr><tr><th rowspan="2">N<sub>s</sub></th><th colspan="3">R<sub>1</sub></th><th colspan="3">R<sub>2</sub></th></tr><tr><th>Opt</th><th>A-EMS</th><th>DMS</th><th>Opt</th><th>A-EMS</th><th>DMS</th></tr></thead><tbody><tr><td>1</td><td>2.70</td><td>2.70</td><td>2.70</td><td>5.15</td><td>5.15</td><td>5.15</td></tr><tr><td>8</td><td>1.45</td><td>1.52</td><td>1.83</td><td>1.85</td><td>2.26</td><td>2.59</td></tr><tr><td>16</td><td>1.32</td><td>1.41</td><td>1.65</td><td>1.52</td><td>1.78</td><td>2.16</td></tr><tr><td>32</td><td>1.19</td><td>1.29</td><td>1.51</td><td>1.29</td><td>1.58</td><td>1.94</td></tr><tr><td>64</td><td>1.12</td><td>1.22</td><td>1.36</td><td>1.08</td><td>1.28</td><td>1.63</td></tr></tbody></table>

### 4.3. Error models investigation

The previous results have shown that A-EMS can select, for most queried values of $q$, a mesh that yields an error comparable to the optimal error and significantly smaller than that obtained with DMS. This selection relies on substituting the true error with the restriction error and predicting values at arbitrary parameters using a GP model. Consequently, the performance of the A-EMS method depends directly on the quality of the GP models.

To assess the validity of this approach, Figure 4a presents a $q-q$ plot comparing exact and restriction errors. Each point corresponds to the restriction error $\|u(q) - \Pi_i u(q)\|$ against the true error $\|u(q) - u_{H_i}(q)\|$, for one of the $N_L$ meshes and one parameter value from the 80 × 80 uniform grid covering $R_2$. The plot thus contains 16 × 80² data points. A strong correlation between exact and restriction errors is evident,
particularly at low error levels, where precision is critical for accurate mesh selection. We also note that while the exact error typically remains above $\epsilon_{ad}$, the restriction error may vanish since it does not incorporate $\epsilon_{ad}$.

Figure 4b compares the restriction errors predicted by the GP models with the computed restriction errors. The correlation remains clear but exhibits greater dispersion than in the previous plot, and the GP models may even predict negative restriction errors. Finally, Figure 4c compares the GP predictions directly with the exact errors. The combined effect of trend errors and the dispersion introduced by restriction and GP modeling does not significantly degrade overall accuracy, thereby explaining the correct behavior of the method.

<figure><img src="image_5.png" /><figcaption>Figure 4: Comparison of exact, restriction, and GP error model. Oblique shock on $R_2$ with $N_L = 16$.</figcaption></figure>

The conclusion from our experiments on the oblique shock problem is that precise error prediction is not required to correctly decide on the mesh from the library to be used at a new query point. Indeed, confusion between meshes yielding comparable errors has little impact on the resulting accuracy. At the same time, even relatively crude GP models can reliably distinguish between well-adapted and poorly adapted meshes.

To further examine the approximation of the error by GP models, Figure 5 analyzes the influence of the neighborhood size considered for GP construction. For simplicity, we adopt a geometric definition of the neighborhood, retaining the closest library parameters $q_j$ to build each GP model $e_i(q)$. The figure reports the dependence of the A-EMS averaged sub-optimality on the neighborhood size for both the large and reduced parametric ranges. The results show that, beyond a sufficient neighborhood size, a minimum sub-optimality is attained. Two distinct behaviors then emerge as the neighborhood size increases. If the library is sufficiently rich to achieve an average error close to $\epsilon_{ad}$, the sub-optimality remains nearly constant and independent of neighborhood size. Conversely, if the library cannot achieve an average error close to $\epsilon_{ad}$, the sub-optimality stays low but increases once the neighborhood size exceeds a certain threshold. This increase is explained by the relatively sparse sampling in low-error regions, which leads to a higher rate of non-optimal selections as the error models become biased by longer-distance trends.

### 4.4. Cost analysis

To further assess the efficiency of A-EMS, we provide a brief analysis of the computational cost and a comparison with other adaptation strategies. Specifically, Figure 6 reports the effective averaged errors as a function of computational cost for the Nominal Mesh Adaptation (NMA), Mean Mesh Adaptation (MMA), A-EMS, and Specific Mesh Adaptation (SMA).

To generate these error-to-cost curves, we progressively decrease the target precision $\epsilon_{ad}$ in the mesh adaptations. We then compute the average error achieved with the adapted meshes and the considered strategies. For SMA, the computational cost is defined by equation (18) as the sum of costs of all adaptation iterations required to reach the average number of iterations $\bar{k}_{ad}$, estimated over the library meshes for each error target. Similarly, the cost of A-EMS is computed using equation (19), corresponding to the cost of
<figure><img src="image_6.png" /><figcaption>Figure 5: Averaged sub-optimality of A-EMS as a function of the neighborhood size in the GP models construction and for the two parametric ranges. Oblique shock on $R_2$ with $N_L = 64$.</figcaption></figure>

solving on the final grid of a $\bar{k}_{ad}$-iteration process. For A-EMS, $\bar{k}_{ad}$ is taken as the weighted average of the selected meshes' $k_{ad}$ according to the selection policy, reflecting potential bias toward specific library meshes. For NMA, the cost corresponds to the last grid of a mesh adapted at the median parameter values, while for MMA the cost derives from solving on a single mesh adapted to $\varepsilon_{ad}$ [13]. In all cases, the same refinement fraction $\alpha_{ad} = 0.08$ is used.

Figure 6 highlights the limitation of the NMA approach: it cannot achieve arbitrary levels of averaged error, because a mesh adapted to a nominal condition is generally unsuitable for other conditions unless the refinement is uniform. Consequently, the NMA cost increases even when the averaged error stagnates. A-EMS exhibits a similar trend, with a rapid initial reduction of the averaged error followed by a plateau, while the cost continues to rise. The asymptotic error level of A-EMS depends on $N_L$ for the same reason as for NMA, which is a particular case of A-EMS with $N_L = 1$. The convergence rate before stagnation is largely independent of $N_L$, as it depends on the average size of the adapted library meshes. In the limit $N_L \to \infty$, A-EMS achieves the same convergence rate as SMA but at a fraction of its computational cost, corresponding to the saved adaptation iterations. In the examples presented, A-EMS is asymptotically five times less expensive than SMA. Nonetheless, very large libraries may introduce storage and handling challenges during mesh selection.

Regarding MMA, we observe that in this example it outperforms SMA at large error levels but does not surpass A-EMS for the considered values of $N_L$. Still, being a single-mesh method, MMA is simpler to implement and still allows for average error control, without encountering convergence plateaus like NMA or A-EMS.

Finally, Figure 7 illustrates the impact of the refinement fraction $\alpha_{ad}$ on A-EMS performance. For $\alpha_{ad} \in [0.02, 0.3]$, the computational cost remains essentially constant for a fixed averaged error. However, excessively large $\alpha_{ad}$ values (e.g., $\alpha_{ad} = 0.6$) reduce adaptation sharpness and therefore the efficiency of A-EMS. Accordingly, it is generally advantageous to use A-EMS with a low refinement fraction.

# 5. Validation on scramjet flow simulation

A two-dimensional supersonic scramjet flow, originally introduced in [16] and widely employed to assess adaptation methods [25, 26, 27, 28], is considered in this section to validate the behavior of EMS.

## 5.1. Scramjet case presentation

The symmetric scramjet domain consists of five regions, as illustrated in Figure 8a: an upstream section, a central duct, two outer ducts (top and bottom), and a downstream section where the duct flows merge. The
<figure><img src="image_7.png" /><figcaption>Figure 6: Cost-efficiency of the different adaptation strategies, using $\alpha_{ad} = 0.08$. Oblique shock on $R_2$.</figcaption></figure>

<caption>Figure 7: Cost-efficiency of the A-EMS method for different refinement fraction. Oblique shock on $R_2$ with $N_L = 16$.</caption>

boundary conditions consist of a supersonic inlet and outlet, with slip conditions imposed on the solid walls. The Euler equations for a perfect gas are solved in the flow domain using the finite-volume, cell-centered solver `elSA` [29] (ONERA-SAFRAN property). Fluxes are discretized using Roe scheme [20], and the unsteady governing equations are time-integrated until steady state using an implicit backward Euler scheme coupled with a LUSSOR-type linear solver.

The variable flow parameters are the inflow Mach number $M_0$ and the inlet incidence angle $\alpha_{in}$, such that $q(\theta) = (M_0(\theta), \alpha_{in}(\theta))$. These two parameters govern the flow structure and determine the position of the shock waves in the domain. The Mach number varies in the range $3 \le M_0 \le 5$, which ensures a constant flow regime and a supersonic outlet across the entire interval. The incidence angle is restricted to $0 \le \alpha_{in} \le 6$, since viscous effects would become non-negligible at higher values. The two parameters are assumed to vary independently, with uniform distributions over their respective ranges.

To apply the EMS method, flow solutions and adapted meshes are computed for $N_L$ flow conditions ($M_0, \alpha_{in}$) generated from a Sobol sequence in the rectangle $[3, 5] \times [0, 6]$. Unless specified otherwise, $N_L = 16$. Each computation is initialized from the same coarse mesh $\mathcal{H}_0$, containing $|\mathcal{H}_0| = 3,000$ cells. Mesh adaptation relies on an interpolation error estimate of the Mach number field to guide refinement, and it terminates when the error estimate $\hat{\varepsilon}_{\mathcal{H}_k}$ falls below $\varepsilon_{ad} = 0.004$. This threshold corresponds approximately to a one-order-of-magnitude reduction of the error estimate compared to the initial mesh $\mathcal{H}_0$.

Figures 8b to 8d show three adapted meshes (upper part of the domain) and their solutions' Schlieren fields (bottom of the domain). The flows present oblique shocks and expansion fans along with complex interactions in the central and outer ducts. The upstream section reduction induces asymmetric oblique shocks that hit the struts past their leading edge. In the outer ducts, their reflections bounce off the walls up to the outlet. The varying duct section creates several expansion fans and oblique shocks interacting with each other in the downstream section. In the central duct, the shocks induced by the struts' leading edges meet and reflect in the converging section. The throat creates a system of interacting expansion fans and shocks that develop and interact in the diverging section down to the outlet. The merging of the central and outer flows and their interaction with the trailing edge structures yield a complex flow in the downstream section. To capture the flow features, the meshes show local adaptations along the shocks and expansion fans in the upstream section and entrance of the central and outer ducts. Further downstream, the diffusion of the expansion fans and the complex interactions lead to more uniform refinement meshes. In addition, comparing the three meshes highlights the displacement and the changes in the flow structures with the conditions.

For the analysis of EMS, a validation set of 100 reference solutions is generated to evaluate the mesh selection strategies and the resulting errors. The validation samples correspond to flow conditions ($M_0, \alpha_{in}$) drawn randomly within their respective ranges. Reference solutions are computed on a fine mesh $\mathcal{H}_{ref}$,
<figure><img src="image_8.png" /><figcaption>Figure 8: Scramjet problem description.</figcaption></figure>

obtained by uniformly refining the initial mesh $\mathcal{H}_0$ six times, resulting in $|\mathcal{H}_{\text{ref}}| \approx 12 \times 10^6$ cells. For consistency, a maximum of six successive refinements is permitted in the AMR procedure, ensuring that all adapted meshes in the EMS library are sub-meshes of $\mathcal{H}_{\text{ref}}$.

## 5.2. Selection maps and error analysis

We first examine the quality of the selection maps for the scramjet problem. The results are presented in Figure 9 for $N_L = 16$. The first plot in Figure 9a illustrates the optimal selection, with segments connecting each validation sample (solid circles) to the adapted mesh in the library (white-core circles) that yields the lowest error in the Mach field. For this case, the error is evaluated using the reference solutions. The plot emphasizes the dominant influence of the Mach number on the optimal selection across all incidence values, which is expected since the Mach number primarily determines the longitudinal position of the shocks, whereas the inlet incidence produces only minor vertical displacements of the shocks within the domain.

Next, we apply the A-EMS method, which relies on Gaussian process models of the restriction errors. Figure 9b displays the discrete map, where green segments indicate selections that match the optimal ones (from the previous figure), and red segments indicate selections that differ from the optimal. Light gray dashed lines indicate the optimal selections for reference. Approximately 80% of the validation samples are associated with their optimal mesh from the library. The A-EMS model successfully captures the dominant influence of the inflow Mach number on the selection process. Moreover, no significant correlation is observed between the selection error and the distance between a sample and its (erroneously) selected mesh, suggesting that the errors in the Gaussian process models exhibit no systematic structure.

The performance of A-EMS contrasts sharply with that of the parametric distance-based selection (DMS), as shown in Figure 9c. Specifically, associating each sample with the library element having the closest flow conditions (accounting for the relative ranges of the parameters) yields the optimal selection for only about 30% of the samples.

In addition to selection performance, the errors achieved by the A-EMS method have been analyzed. Figure 10 shows histograms of the true errors for the optimal selection (black), A-EMS (blue), and DMS (red), along with the corresponding estimated distributions. We first observe that the error density of A-EMS closely matches that of the optimal selection, albeit with a slightly heavier tail toward higher error values. This is expected, as A-EMS correctly selects the optimal mesh for the majority of samples. The comparison of the densities indicates, however, that non-optimal selections do not lead to a substantial increase in
<figure><img src="image_9.png" /><figcaption>Figure 9: Selection maps for the scramjet problem with $N_L = 16$. The maps show the association between the validation samples (full circles) and an element of the library (white-core circles). The optimal associations are shown with green segments, and non-optimal are shown in red.</figcaption></figure>

the true error; these miss-selected meshes remain effective choices, though not optimal. Quantitatively, A-EMS exhibits a sub-optimality of less than 1% relative to the optimal selection. In contrast, the true error distribution for the distance-based selection is considerably flatter and exhibits a pronounced high-error tail, corresponding to a sub-optimality of 26.5%. Increasing the library size to $N_L = 32$ reduces the optimal averaged error from 0.0061 to 0.0053; A-EMS again approaches this optimal error within 1% sub-optimality, while the distance-based method improves but remains 24% sub-optimal.

<figure><img src="image_10.png" /><figcaption>Figure 10: True error distribution and histogram using A-EMS, DMS, and optimal selections. Scramjet case with $N_L = 16$.</figcaption></figure>

The error analysis above indicates that the averaged error ($\approx 0.0061$) exceeds the target error $\varepsilon_{ad} = 0.004$ for all selection approaches, suggesting that $N_L = 16$ is insufficient to fully cover the parameter domain.

To complete the analysis, we compare with alternative approaches. A naive single-mesh approach based on a hierarchical mesh containing all library meshes (the union mesh) would require roughly $10^6$ cells, which is computationally prohibitive. Alternatively, the MMA method [13] can generate a single mesh capable of achieving an averaged error comparable to the optimal selection (0.0061), with $|\mathcal{H}_{MMA}| = 3 \times 10^5$ cells approximately, as size 50% larger than the average size of the adapted meshes in the library. Furthermore, to ensure an average error close to the adaptation threshold $\varepsilon_{ad} = 0.004$, MMA requires a mesh with $5 \times 10^5$ cells, an even bigger overshoot of the library’s average mesh size. In contrast, by increasing the library density, A-EMS can reach $\varepsilon_{ad}$ without increasing the average mesh size. Therefore, for this problem and the
considered error level, the A-EMS method enables the use of meshes that are at least 2–3 times coarser than a single-mesh approach, provided that the library densely samples the parameter domain.

### 5.3. Error models investigation

To better assess the effectiveness of A-EMS, Figure 11 presents an analysis of the predicted errors. The plot shows the GP-predicted errors at the validation samples as a function of the true error for each mesh in the library, yielding $N_L \times 100 = 1,600$ data points. As expected, the predicted errors tend to underestimate the true errors, since they are constructed from the restriction errors and neglect additional contributions (see equation (10)). In particular, the GP models predict errors below $\epsilon_{ad}$, which is impossible for the true errors. Two distinct trends are observed depending on the true error magnitude. For errors smaller than $\approx 0.015$, there is a reasonable correlation between the GP predictions and the true errors. This correlation deteriorates for larger errors, and the predictive capability of the GP models quickly degrades. Nevertheless, A-EMS selects sufficiently close meshes to ensure that the true error rarely exceeds 0.01, as seen in Figure 10. Thus, the GP models provide rough error estimates that are adequate to correctly guide the mesh selection from the library.

<figure><img src="image_11.png" /><figcaption>Figure 11: Predicted error of A-EMS (GP models) against true error. Scramjet case with $N_L = 16$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 12: Evolution of the sub-optimality measure with the size of the neighborhood used for A-EMS GP models construction. Scramjet case with $N_L = 32$.</figcaption></figure>

Finally, Figure 12 illustrates the dependence of the sub-optimality $\bar{\rho}_{A-EMS}$ on the neighborhood size used to construct the Gaussian process models. Here, $N_L = 32$ is employed to allow for a sufficient range of neighborhood sizes. The observed trends are consistent with those reported for the oblique shock case. For small neighborhoods, the sub-optimality is high but decreases rapidly as the neighborhood size increases. Beyond a certain size, the sub-optimality stagnates and then gradually rises. This behavior reflects the trade-off between GP models designed to capture local error structures versus those optimized for a more global fit. These results highlight the importance of carefully choosing the neighborhood size to balance local accuracy with the broader generalization capability of the GP models.

## 6. Conclusion

We have introduced the Error-based Mesh Selection method, A-EMS, to improve the efficiency of simulations performed across varying conditions. The method predicts, for any given condition, the discretization error that would result from solving on each mesh in a precomputed library, and selects the mesh expected to minimize this error. These predictions rely on Gaussian process models of restriction errors, enabling error estimation across all library meshes without additional model solves. The approach was assessed on
an analytical oblique shock problem and on a finite-volume CFD configuration involving a supersonic flow. In both cases, the restriction error proved to be a reliable surrogate for the discretization error, and the Gaussian process models enabled mesh selections close to the optimal choice, significantly outperforming simpler strategies such as parametric distance-based selection. While achieving low error levels requires a sufficiently large library, the complexity analysis confirms that this requirement does not undermine the overall efficiency of the approach. A-EMS consistently outperforms nominal mesh strategies and is more efficient than repeated condition-specific adaptations, provided the library contains enough adapted meshes to meet the desired accuracy.

Despite these advantages, several challenges remain. Gaussian process modeling becomes difficult when the restriction error exhibits strong anisotropy or when the library is unevenly distributed across the parameter space. Sequential construction of the library is, therefore, a promising direction. Future work will focus on designing infilling criteria guided by GP uncertainty, to quantify mis-selection likelihood and suboptimality risk. More advanced GP modeling techniques may further improve the selection quality. Additionally, coupling A-EMS with Mean Mesh Adaptation (MMA) to compute a mean adapted mesh within each selection region represents a potential avenue to enhance robustness.

## Acknowledgment

The authors would like to thank Niels Brun for his contribution to the development of the library creation and assessment tools used in this work.

This work was partially supported by the SONICE project, which is a France Relance project co-funded by the French Civil Aviation Authority (DGAC). France Relance benefits from EU funding via the NextGenerationEU initiative.

## References

[1] F. Alauzet, A. Loseille, A decade of progress on anisotropic mesh adaptation for computational fluid dynamics, Computer-Aided Design 72 (2016) 13-39.

[2] W. G. Habashi, J. Dompierre, Y. Bourgault, D. Ait-Ali-Yahia, M. Fortin, M.-G. Vallet, Anisotropic mesh adaptation: Towards user-independent, mesh-independent and solver-independent CFD. Part I: General principles, International Journal for Numerical Methods in Fluids 32 (2000) 725-744. doi:doi: 10.1002/(SICI)1097-0363(20000330)32:6<725::AID-FLD935>3.0.CO;2-4.

[3] T. J. Baker, Mesh adaptation strategies for problems in fluid dynamics, Finite Elements in Analysis and Design 25 (1997) 243-273.

[4] F. Alauzet, L. Frazza, D. Papadogiannis, Periodic adjoints and anisotropic mesh adaptation in rotating frame for high-fidelity RANS turbomachinery applications, Journal of Computational Physics 450 (2022) 110814.

[5] G. Vivarelli, N. Qin, S. Shahpar, D. Radford, Efficient adjoint-based mesh adaptation applied to turbo-machinery flows, in: Proceedings of ASME Turbo Expo 2018, 2018, pp. 1-14.

[6] H. Bijl, Uncertainty Quantification in Computational Fluid Dynamics, number v.92 in Lecture Notes in Computational Science and Engineering Ser, 1st ed., Springer International Publishing AG, Cham, 2013.

[7] J. M. Burt, M. A. Hagenmaier, D. R. Eklund, H. L. Wilkin, E. Josyula, Uncertainty Quantification for a Scramjet Inlet Flow, in: 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, American Institute of Aeronautics and Astronautics, Cleveland, OH, 2014, pp. 1-13. doi:doi: 10.2514/6.2014-3844.

[8] B. Zhang, Y. Ju, Z. Li, X. Zhang, Y. Liu, C. Zhang, Reconstruction and reduction of realistic manufacturing error field for uncertainty quantification of transonic axial compressor rotor, Journal of Turbomachinery 148 (2026) 11002. doi:doi: 10.1115/1.4069132.

[9] J. Van Langenhove, D. Lucor, F. Alauzet, A. Belme, Goal-oriented error control of stochastic system approximations using metric-based anisotropic adaptations, Journal of Computational Physics 374 (2018) 384-412. doi:doi: 10.1016/j.jcp.2018.07.044.

[10] M. Capriati, A. Cortesi, T. E. Magin, P. M. Congedo, Stagnation point heat flux characterization under numerical error and boundary conditions uncertainty, European Journal of Mechanics / B Fluids 95 (2022) 221-230.

[11] F. Palacios, K. Duraisamy, J. J. Alonso, E. Zuazua, Robust Grid Adaptation for Efficient Uncertainty Quantification, AIAA Journal 50 (2012) 1538-1546. doi:doi: 10.2514/1.J051379.

[12] N. Barral, T. Taddei, I. Tifouti, Registration-based model reduction of parameterized PDEs with spatio-parameter adaptivity, Journal of Computational Physics 499 (2024) 112727. doi:doi: 10.1016/j.jcp.2023.112727.

[13] H. Dornier, O. Le Maître, P. Congedo, I. Salah El Din, J. Marty, S. Bourasseau, Mean Mesh Adaptation for Efficient Cfd Simulations with Operating Conditions Variability, 2024. doi:doi: 10.2139/ssrn.5039807.
[14] A. Gersho, R. M. Gray, Vector Quantization and Signal Compression, Springer US, Boston, MA, 1992. doi: 10.1007/978-1-4615-3626-0.

[15] R. Gray, Vector quantization, IEEE ASSP Magazine 1 (1984) 4–29. doi:doi: 10.1109/MASSP.1984.1162229.

[16] A. Kumar, Numerical Analysis of the Scramjet-Inlet Flow Field by Using Two-Dimensional Navier-Stokes Equations, Technical Report NASA-TP-1940, NASA Langley Research Center, 1981.

[17] R. Löhner, K. Morgan, O. Zienkiewicz, An adaptive finite element procedure for compressible high speed flows, Computer Methods in Applied Mechanics and Engineering 51 (1985) 441–465. doi:doi: 10.1016/0045-7825(85)90042-8.

[18] J. T. Oden, T. Strouboulis, Ph. Devloo, Adaptive finite element methods for high-speed compressible flows, International Journal for Numerical Methods in Fluids 7 (1987) 1211–1228. doi:doi: 10.1002/fld.1650071105.

[19] J. Peraire, M. Vahdati, K. Morgan, O. Zienkiewicz, Adaptive remeshing for compressible flow computations, Journal of Computational Physics 72 (1987) 449–466. doi:doi: 10.1016/0021-9991(87)90093-3.

[20] P. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics 43 (1981) 357–372. doi:doi: 10.1016/0021-9991(81)90128-5.

[21] P. Wesseling, An Introduction to Multigrid Methods, Pure and Applied Mathematics, Wiley, Chichester, 1995.

[22] C. E. Rasmussen, C. K. I. Williams, Gaussian Processes for Machine Learning, Adaptive Computation and Machine Learning, 3. print ed., MIT Press, Cambridge, Mass., 2008.

[23] P. Saves, R. Lafage, N. Bartoli, Y. Diouane, J. Bussemaker, T. Lefebvre, J. T. Hwang, J. Morlier, J. R. Martins, SMT 2.0: A Surrogate Modeling Toolbox with a focus on hierarchical and mixed variables Gaussian processes, Advances in Engineering Software 188 (2024) 103571. doi:doi: 10.1016/j.advengsoft.2023.103571.

[24] J. D. Anderson, Fundamentals of Aerodynamics, McGraw-Hill Series in Aeronautical and Aerospace Engineering, sixth edition ed., McGraw-Hill Education, New York, NY, 2017.

[25] M. J. Castro-Díaz, F. Hecht, B. Mohammadi, O. Pironneau, Anisotropic unstructured mesh adaption for flow simulations, International Journal for Numerical Methods in Fluids 25 (1997) 475–491. doi:doi: 10.1002/(SICI)1097-0363(19970830)25:4<475::AID-FLD575>3.0.CO;2-6.

[26] F. Alauset, Adaptation de Maillage Anisotrope En Trois Dimensions. Application Aux Simulations Instationnaires En Mécanique Des Fluides., Ph.D. thesis, Montpellier II, 2003.

[27] A. Loseille, A. Dervieux, P. Frey, F. Alauzet, Achievement of global second order mesh convergence for discontinuous flows with adapted unstructured meshes, American Institute of Aeronautics and Astronautics (2007).

[28] J. Van Langenhove, Adaptive Control of Deterministic and Stochastic Approximation Errors in Simulations of Compressible Flow, Ph.D. thesis, Université Pierre et Marie Curie - Paris VI, Paris, 2017.

[29] L. Cambier, S. Heib, S. Plot, The ONERA elsA CFD software : Input from research and feedback from industry, Mechanics & Industry, EDP Sciences 14 (2013) 159–174.