# Enhancing physical consistency in stochastic optimization for adjoint-based inverse problems: Application to compressible RANS simulations in the discontinuous Galerkin framework

Bartolomeo Fanizza, Pedro Stefanin Volpiani, Florent Renac, Denis Sipp

► To cite this version:

Bartolomeo Fanizza, Pedro Stefanin Volpiani, Florent Renac, Denis Sipp. Enhancing physical consistency in stochastic optimization for adjoint-based inverse problems: Application to compressible RANS simulations in the discontinuous Galerkin framework. Applied Mathematical Modelling, 2026, 153, pp.116606. 10.1016/j.apm.2025.116606. hal-05465442

HAL Id: hal-05465442

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

Submitted on 19 Jan 2026

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 CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Enhancing physical consistency in stochastic optimization for adjoint-based inverse problems: Application to compressible RANS simulations in the discontinuous Galerkin framework

Bartolomeo Fanizza ID a,b,*, Pedro Stefanin Volpiani ID c, Florent Renac ID a,
Denis Sipp ID c

## ARTICLE INFO

### Keywords:

- Stochastic optimization
- ADAM
- BFGS
- Data assimilation
- Inverse problems
- RANS equations
- Turbulence
- Discontinuous Galerkin

## ABSTRACT

This work explores the use of quasi-Newton and stochastic optimization methods for gradient-based inverse problems in physical modeling, focusing on their application within high-order numerical frameworks. Such problems are often characterized by the following challenges: (i) physical constraints often yield a highly non-convex design space with multiple local optima; (ii) solutions must satisfy intrinsic properties, such as boundary and regularity conditions, which are not easily enforced as explicit constraints. Conventional stochastic optimizers, such as N-ADAM, exhibit significant information loss in their update rules, leading to non-physical solutions. To overcome these limitations, we introduce a novel stochastic optimizer, V-N-ADAM-DG, which incorporates adjoint-state information into the update rule to maintain physically meaningful corrections in terms of regularity and boundary conditions. We validate our approach in the context of mean-flow reconstruction for Reynolds-averaged Navier-Stokes (RANS) simulations using a high-order discontinuous Galerkin (DG) discretization method, as proposed by Fanizza et al. (2025). The optimization framework considers both vectorial corrective terms, inferred in the momentum and energy equations, and scalar corrective terms in the Spalart-Allmaras (SA) transport equation. The V-N-ADAM-DG optimizer effectively reconstructs mean flow quantities while ensuring smooth transitions of corrective parameters at boundaries, an improvement over standard stochastic optimizers. Additionally, it facilitates a rapid decay of the optimal degrees of freedom (DOFs), leading to smoother corrections in high-order reconstructions-achieving a balance between the robustness of quasi-Newton methods (such as L-BFGS) and the flexibility of stochastic approaches. Numerical experiments across various flow configurations demonstrate that V-N-ADAM-DG consistently outperforms both L-BFGS and N-ADAM, particularly in complex inverse problems that employ multiple combined cost functions to reconstruct different physical quantities.

## 1. Introduction

Computational optimization is an interdisciplinary field at the intersection of computer science and applied mathematics. It involves the use of algorithms and computational methods to find the optimal solution to a given problem by selecting a set of

<footnote>a Department of Aerodynamics, Aeroelasticity and Acoustics (DAAA), ONERA, 29 Av. de la Division Leclerc, Châtillon, 92320, France</footnote>

<footnote>b Ecole Polytechnique, Institut Polytechnique de Paris, Route de Saclay, Palaiseau, 91120, France</footnote>

<footnote>c Department of Aerodynamics, Aeroelasticity and Acoustics (DAAA), ONERA, 8 Rue des Vertugadins, Meudon, 92190, France</footnote>

<footnote>* Corresponding author.
E-mail address: bartolomeo.fanizza@onera.fr (B. Fanizza).</footnote>
parameters or decisions from a range of possibilities while minimizing a predefined cost function. This field is widely applied across various disciplines, including engineering and industry, ranging from structural mechanics [1] to aeronautical engineering design [2,3].

Computational optimization is also employed to improve the accuracy of low-fidelity models, often expressed in the form of partial differential equations (PDEs), by incorporating information from high-fidelity solutions. By using measurements from high-fidelity solutions, an optimization problem can be formulated to minimize the discrepancy between low- and high-fidelity observed quantities through the tuning of a correction term. This correction term is inferred in the data-augmented (DA)-PDE problem, which must be satisfied. The DA-PDE problem serves as the constraint of the optimization problem, often referred to as an *inverse problem*.

Once the constrained optimization problem is defined, search algorithms are employed to determine the optimal solution. A widely used class of search algorithms consists of gradient-based optimization methods. These algorithms utilize information about the objective function and its gradient with respect to the tuning parameters to infer the optimal correction for the low-fidelity model. When the set of corrective parameters is large, an adjoint formulation can be introduced to efficiently compute the required gradients. The combination of these search algorithms with the adjoint formulation constitutes a powerful approach for solving inverse problems. In computational fluid dynamics (CFD), optimization methods increasingly rely on adjoint-based techniques to efficiently compute gradients in high-dimensional parameter spaces. Jameson (1988) [4] demonstrated the effectiveness of control-theory and adjoint formulations for aerodynamic design, establishing the foundation for gradient-based flow optimization. Later, Giles and Pierce (2000) [5] introduced the adjoint approach to design optimization and sensitivity analysis in CFD. More recently, Martins and Hwang (2013) [6] reviewed derivative computation methods for large-scale multidisciplinary models, highlighting the ongoing central role of adjoint methods in modern engineering optimization, and Li and Zheng (2017) [7] summarized how turbo-machinery blade-design workflows employ stochastic algorithms to escape local minima.

A relevant application, in line with the topic of this work, involves the use of optimization techniques to refine Reynolds-averaged Navier-Stokes (RANS) simulations through field inversion (FI) [8,9]. Due to their low computational cost, RANS simulations are widely preferred in aerodynamic system design. However, RANS models exhibit significant shortcomings when applied to flows with separation, high-streamline curvature, strong pressure gradients, or turbulent mixing.

Several alternative corrective formulations have been proposed. Duraisamy et al. [8] formulated a FI problem, leveraging high-fidelity measurements to improve RANS simulations assisted by finite-volumes methods. They optimized a corrective parameter defined in physical space to tune the rate of turbulence production. Franceschini et al. [10] optimized source terms in the momentum equations and in the Spalart-Allmaras model to enhance RANS predictions for strongly detached incompressible flows. Cato et al. [11] compared six different DA techniques across various flow configurations. They improved the source term in the Spalart-Allmaras model by replacing it with a modified source term scaled by the turbulent variable, ensuring that corrections were concentrated in regions of high turbulence intensity. Ferrero et al. [12] and Fanizza et al. [13] investigated field-inversion problems respectively for transonic and supersonic flows using a high-order discretization method based on the discontinuous Galerkin (DG) formalism. In [13], various corrective formulations, discretized in the DG space, have been employed to reconstruct the mean flow in a shock-wave/boundary layer interaction (SWTBLI) configuration featuring a strong recirculation bubble. The discrete adjoint problem was formulated to compute the gradient during optimization. This approach successfully reconstructed the mean velocity field and wall quantities, such as pressure and skin friction coefficients, by introducing corrections in the momentum and energy equations that targeted the modeled divergence of the Reynolds stress tensor (RST) and turbulent heat flux. Notably, all previous works have relied on the L-BFGS algorithm or other deterministic optimizers to solve the optimization problem. L-BFGS is a widely used quasi-Newton method [14,15] that estimates the inverse Hessian matrix to guide its search through variable space, making it well-suited for large-scale problems. However, while quasi-Newton methods are highly effective as local search algorithms, they may struggle to find the global optimum if the problem is non-convex. An alternative class of optimization algorithms consists of stochastic gradient descent methods. The ADAM optimizer [16] is a prominent example within this class. Based exclusively on first-order gradient information, it has low memory requirements and computes adaptive learning rates for each parameter using the first and second moments of the gradients. Although stochastic optimizers have achieved considerable success in machine learning and deep learning—particularly in training artificial neural networks (ANNs)—their application in FI-based approaches remains limited.

This work focuses on the use of stochastic optimizers to solve adjoint-based complex inverse problems within the DG numerical framework or similar high-order numerical methods. The optimal degrees of freedom (DOFs), used for the numerical discretization of the corrective parameters, are sought in the DG space, with a discrete adjoint approach employed to compute the gradient of the cost function in a manner consistent with the formal problem. The identification of specific properties of the correction term that ensure accurate boundary and asymptotic behaviors of the modeled terms is directly linked to the adjoint solution. These properties are analyzed within the context of the optimizer's update rule. While quasi-Newton methods naturally preserve these properties, stochastic algorithms introduce additional challenges. To address these challenges, a key contribution of this work is the development of a new stochastic search algorithm designed to manage highly non-convex objective functions while embedding physical constraints into the optimal corrective term through a more stable update rule. The modifications are inspired by the vectorial version of the ADAM algorithm proposed by Ling et al. [17] for geometrical optimization problems. Here, additional refinements are introduced to enhance the smoothness of the corrective parameters represented in the DG space, ensuring a rapid decay of modal coefficients in high-order polynomial representations. We demonstrated the benefits of the proposed modifications in the application of mean flow reconstruction for RANS simulations, assisted by a modal DG method. We considered both vectorial corrective terms, inferred in the momentum and energy equations, and scalar corrective terms in the Spalart-Allmaras (SA) transport equation. The FI problem was solved using combined cost functions, which highlighted the limitations of the L-BFGS and N-ADAM optimizers. The novel stochastic
optimizer successfully achieved lower minima than the L-BFGS optimizer when handling complex cost functions, while also ensuring feasible optimal solutions that preserved the physical properties of the RANS equations and the turbulence model.

The paper is structured as follows: Section 2 introduces the inverse-problem methodology for data-augmented numerical simulations, including the discrete adjoint formulation and gradient-based optimization. Section 3 examines the inferred correction's physical role and analyzes optimizer update rules for quasi-Newton and stochastic optimizers are analyzed to identify potential shortcomings in the physical framework of the inverse problem. Section 4 applies the methodology to mean-flow reconstruction in RANS simulations, introducing two formulations for data assimilation. Section 5 evaluates the performance of the optimizer in inverse problems of laminar flows using high-order polynomials in the DG representation, as well as in inverse problems aimed at reconstructing various mean-flow quantities in the SWTBLI configuration. Finally, Section 6 summarizes key conclusions.

## 2. Inverse problems for data-augmented physical models

### 2.1. Baseline and data-augmented problems

Consider the generic physical model described by the evolution of the solution $u(x, t)$, in a space domain $\Omega \subset \mathbb{R}^d$ and time, governed by $N_{eq}$ nonlinear equations, often in the form of PDEs:

$$ \mathcal{R}_h(u(x, t)) = 0 \quad \text{in} \quad \Omega \times (0, +\infty) \qquad (1) $$

with appropriate initial and boundary conditions on $\partial\Omega$. To solve this problem, the computational domain $\Omega$ is discretized into a shape-regular mesh $\Omega_h$, consisting of $N_{elts}$ non-overlapping elements $\kappa$. The numerical solution $u_h$ is then sought in a specific function space, determined by the chosen numerical scheme, using the so-called degrees of freedom (DOFs) U. The steady-state problem (1) is rewritten in discrete form as

$$ \mathbf{R}(U) = 0. \qquad (2) $$

The numerical solution $u_h$ is then reconstructed from the solution DOFs of (2). The accuracy of the approximate solution may be improved via a functional correction term $f_{DA}(x)$ defined in the physical space and introduced in problem (1). This correction can be optimized using available measurements $m$ of the ground-truth or high-fidelity solution $u_m$. The functional correction must also be discretized within the computational mesh $\Omega_h$ through the DOFs vector $F_{DA}$. Consequently, the so-called data-augmented (DA) problem is expressed as

$$ \mathbf{R}(U, F_{DA}) = 0. \qquad (3) $$

### 2.2. Numerical discretization

For the sake of illustration, we consider a discontinuous Galerkin (DG) formulation to solve both the baseline and DA problems. However, the investigations in this work can be easily extended to a wider range of numerical approaches, such as finite-volume methods.

The DG method constructs a discrete weak formulation of the governing equations. The numerical solution is sought within the function space of discontinuous polynomials, defined as

$$ \mathcal{V}_h^p = \{ \phi \in L^2(\Omega_h) : \phi|_{\kappa} \circ F_{\kappa} \in \mathcal{P}^p(I^d), \forall \kappa \in \Omega_h \}, \qquad (4) $$

where $\mathcal{P}^p(I^d)$ denotes the space of functions constructed from tensor products of polynomials of total degree at most $p$ on the reference element $I^d$, with $I = [-1, 1]$. Each physical element $\kappa$ is obtained as the image of $I^d$ under the mapping $F_{\kappa}$. The numerical approximation of the PDE solution is then given by

$$ \mathbf{u}_h(x) = \sum_{l=1}^{N_p} \phi_{\kappa}^{l}(x) U_{\kappa}^{l}, \quad \forall \mathbf{x} \in \kappa, \quad \kappa \in \Omega_h, \qquad (5) $$

where $(U_{\kappa}^l)_{1 \le l \le N_p}$ represents the degrees of freedom (DOFs) within element $\kappa$. The set $(\phi_{\kappa}^l, \dots, \phi_{\kappa}^{N_p})$ forms a hierarchical and orthogonal modal basis of $\mathcal{V}_h^p$, restricted to $\kappa$. The dimension of this space is given by $N_p = \frac{(p+d)!}{p!d!}$.

An orthonormal basis is constructed using a modified Gram-Schmidt orthonormalization procedure, ensuring a diagonal mass matrix $M = \text{diag}(|\kappa| : \kappa \in \Omega_h)$. Further details on this approach can be found in [18,19].

Similarly, the inferred correction term can be represented within the function space of discontinuous polynomials as

$$ \mathbf{f}_{DA}(x) = \sum_{l=1}^{N_p} \phi_{\kappa}^{l}(x) \mathbf{F}_{DA,\kappa}^{l}, \quad \forall \mathbf{x} \in \kappa, \quad \kappa \in \Omega_h, \qquad (6) $$

where the same polynomial order $p$ is used for both the numerical solution and the corrective parameter. The coefficients $(\mathbf{F}_{DA,\kappa}^l)_{1 \le l \le N_p}$ represent the DOFs for the DG representation of the correction term $f_{DA}$.
## 2.3. Inverse problem, adjoint state & gradient computation

The computation of the correction DOFs $\mathbf{F}_{\text{DA}}$ is performed by solving an inverse problem, formulated as a constrained optimization problem:

$$ \begin{array}{l} \mathbf{F}_{\text{DA}}^{\text{opt}} = \underset{\mathbf{F}_{\text{DA}} \in \mathbb{R}^{N_{\text{DOFs}} \times N_r}}{\arg \min} \mathcal{J}_h(\mathbf{U}, \mathbf{F}_{\text{DA}}) \\ \text{s.t.} \quad \mathbf{R}(\mathbf{U}, \mathbf{F}_{\text{DA}}) = 0 \end{array} \qquad (7) $$

where $\mathcal{J}_h(\mathbf{U}, \mathbf{F}_{\text{DA}})$ is the cost function to be minimized, and $\mathbf{R}(\mathbf{U}, \mathbf{F}_{\text{DA}}) = 0$ represents the constraint. The size of the optimization problem is $N_{\text{DOFs}}^{\text{DA}} = N_{\text{DOFs}} \times N_r$, where $N_{\text{DOFs}} = N_{\text{elts}} \times N_p$ denotes the total number of DOFs for the DG representation of one of the $N_r$ components of the correction vector $\mathbf{f}_{\text{DA}}$.

A general cost function would consist of two terms: the measurement error $\mathcal{M}$, which quantifies the discrepancy between the computed solution $\mathbf{u}_h$ and reference data $\mathbf{m}$ in the measurement space $\Omega_m$, and a regularization term $\mathcal{P}$, which may depend on the correction field in the regularization space $\Omega_{\text{reg}}$:

$$ \mathcal{J}_h(\mathbf{U}, \mathbf{F}_{\text{DA}}) = \frac{1}{2} \|\mathcal{M}(\mathbf{u}_h, \mathbf{m})\|_{L^2(\Omega_m)}^2 + \frac{\alpha_p}{2} \|\mathcal{P}(\mathbf{f}_{\text{DA}})\|_{L^2(\Omega_{\text{reg}})}^2 \qquad (8) $$

where we assume the regularization satisfies $\mathcal{P}(0) = 0$ to have norm-like properties. The coefficient $\alpha_p$ regulates the intensity of the penalization term.

Multiple error metrics $\mathcal{M}(\mathbf{u}_h, \mathbf{m})$ are considered in this work. Their expressions depend on the specific physical model being augmented with reference data and are presented in Section 5 alongside the numerical results. For completeness, the linearization of these metrics is discussed in [13].

We combine the high-order DG solver and a general gradient-based optimization algorithm. The gradient of the objective function with respect to the parameters is then required and finite-difference techniques are highly expensive for our test cases. Hence, following the work of Franceschini et al. [20], we adopt the so-called *adjoint approach*, based on the Lagrangian formalism, to compute the required gradient. This formalism was later extended to the DG framework in [13], and the main steps are recalled here for clarity. Further studies on the adjoint consistency analysis of DG discretization are addressed in [21,22].

We introduce the Lagrangian multipliers $\mathbf{z}_h = \sum_{l=1}^{N_p} \Phi_{\mathbf{k}}^{l}(\mathbf{x}) \mathbf{Z}_{\mathbf{k}}^{l}$, in $(\mathcal{V}_h^p)^{d+3}$, associated to the discrete steady DA problem (3) and consider the following Lagrangian functional

$$ \mathcal{L}_h(\mathbf{U}, \mathbf{Z}, \mathbf{F}_{\text{DA}}) = \mathcal{J}_h(\mathbf{U}, \mathbf{F}_{\text{DA}}) - \mathbf{Z}^{\top} \cdot \mathbf{R}(\mathbf{U}, \mathbf{F}_{\text{DA}}), \qquad (9) $$

with $\mathbf{Z}$ the vector of DOFs in $\mathbf{z}_h$. It is possible to show that if the triplet $(\mathbf{U}, \mathbf{Z}, \mathbf{F}_{\text{DA}})$ is such that $\frac{\partial \mathcal{L}_h}{\partial \mathbf{Z}^\xi} = \frac{\partial \mathcal{L}_h}{\partial \mathbf{U}_\mathbf{k}^\xi} = 0$, for $\mathbf{k} \in \Omega_h$ and $1 \le \xi \le N_p$, then

$$ \mathbf{R}(\mathbf{U}, \mathbf{F}_{\text{DA}}) = 0 \qquad (10) $$

$$ \frac{\partial \mathbf{R}^{\top}}{\partial \mathbf{U}_{\mathbf{k}}^{\xi}} \cdot \mathbf{Z} = \frac{\partial \mathcal{J}_{\mathbf{h}}}{\partial \mathbf{U}_{\mathbf{k}}^{\xi}} \qquad (11) $$

Eq. (10) imposes the constraint of the optimization problem (7) to be satisfied. Eq. (11) defines the adjoint state $\mathbf{Z}$ through a system of linear equations involving the Jacobian matrix and a right-hand side linked to the cost function.

Once $\mathbf{U}$ and $\mathbf{Z}$ are computed, the increment of the cost-functional $\delta\mathcal{J}_h$ may be obtained through:

$$ \delta\mathcal{J}_h = \frac{\partial \mathcal{L}_h}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}} \delta\mathbf{F}_{\text{DA},\kappa}^{\xi} = \left\{ \mathbf{M}^{-1} \left[ \frac{\partial \mathcal{J}_h}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}} - \frac{\partial \mathbf{R}^{\top}}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}} \cdot \mathbf{Z} \right] \mathbf{M} \delta\mathbf{F}_{\text{DA},\kappa}^{\xi} \right\}^{\top} \qquad (12) $$

where $\mathbf{M}$ is the mass matrix that takes into account the volume of the elements. The gradient $\nabla_{\mathbf{F}_{\text{DA},\kappa}^{\xi}} \mathcal{J}_h$ inside the above parenthesis is based on the continuous inner-product and should converge to the continuous gradient as the discretization is refined. Yet, in view of using this gradient in search-algorithm routines, which optimizes a cost-functional with the aid of gradients based on canonical inner-products, we rewrite the increment of the cost-functional with a change of variable:

$$ \delta \mathbf{F}_{\text{DA},\kappa}^{\xi} = \mathbf{L}^{-1} \delta \tilde{\mathbf{F}}_{\text{DA},\kappa}^{\xi} \qquad (13) $$

where we have used the Cholesky decomposition of the mass matrix: $\mathbf{M} = \mathbf{L}^{\top}\mathbf{L}$. This decomposition is straightforward here as the mass matrix is diagonal (see Section 2.2). We then obtain:

$$ \delta\mathcal{J}_h = \left( \nabla_{\mathbf{F}_{\text{DA},\kappa}^{\xi}} \mathcal{J}_h \right)^{\top} \delta\tilde{\mathbf{F}}_{\text{DA},\kappa}^{\xi}, \qquad (14) $$
where:

$$ \nabla_{\tilde{\mathbf{F}}_{\text{DA},\kappa}^{\xi}} J_h = (\mathbf{L}^{\top})^{-1} \left[ \frac{\partial J_h}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}} - \frac{\partial \mathbf{R}}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}}^{\top} \cdot \mathbf{Z} \right]. \quad (15) $$

Hence it is the gradient (15) that will be provided to the optimizer, which works internally with the variables $\tilde{\mathbf{F}}_{\text{DA},\kappa}^{\xi}$. To come back to the initial control parameter $\mathbf{F}_{\text{DA},\kappa}^{\xi}$, the output from the optimizers need to be transformed back with Eq. (13). For further details over the linearization of different terms in previous systems we refer to [13].

# 3. Stochastic optimization for inverse problems and physical consistency

## 3.1. Physical significance of the inverse problem solution

The optimal solution to the inverse problem (7) is obtained using iterative search algorithms. Accordingly, the solution can be expressed as an initial guess augmented by $N_{\text{iter}}$ iterative corrections:

$$ \mathbf{F}_{\text{DA}}^{\text{opt}} = \mathbf{F}_{\text{DA}}^{(t=0)} + \sum_{t=1}^{N_{\text{iter}}} \delta \mathbf{F}_{\text{DA}}^{(t)}, $$

where each correction $\delta \mathbf{F}_{\text{DA}}^{(t)}$ at iteration $t$ depends on the optimization algorithm employed and the form of its update rule. In this study, we use updates based on the gradient $g^{(t)} = \nabla_{\tilde{\mathbf{F}}_{\text{DA},\kappa}^{\xi}} J_h$, which naturally inherits the regularity constraints and boundary conditions imposed by the governing equations. Specifically, the correction is defined as a function of the gradient:

$$ \delta \mathbf{F}_{\text{DA}}^{(t)} := \mathbf{F}_{\text{DA}}^{(t)} - \mathbf{F}_{\text{DA}}^{(t-1)} = \mathcal{F}(g^{(t)}). $$

The gradient $g^{(t)}$, as defined in Eq. (15), consists of two main components. First, the transpose of the linearized residual with respect to the tuning parameter, $(\frac{\partial \mathbf{R}}{\partial \mathbf{F}_{\text{DA},\kappa}^{\xi}})^{\top}$, which depends on the structure of the inferred correction in the DA equations. Second, it involves the adjoint state $\mathbf{Z}$, which quantifies the sensitivity of the objective function to perturbations in the tuning variable, and is therefore influenced by the observed measurements.

It is crucial to choose a search algorithm whose update rule $\mathcal{F}$ ensures a smooth dependence on the gradient. This preserves as much gradient information as possible throughout the optimization process, resulting in an optimal field that satisfies physical constraints which would otherwise be difficult to enforce.

Here, we illustrate through an example the significance of the inverse problem solution and the optimized corrective term. We consider the Reynolds-averaged Navier-Stokes (RANS) equations together with a vectorial correction term $\mathbf{f}_{xyz}$ introduced into the momentum-conservation equation, thereby obtaining the so-called DA-RANS equations

$$ \nabla \cdot (\rho \mathbf{v}^{\top} \mathbf{v}) = -\nabla \mathbf{p} + \nabla \cdot (\boldsymbol{\tau} + \boldsymbol{\tau}_t) + \rho \mathbf{f}_{xyz} \quad (16) $$

where $\rho$ is the density, $\mathbf{v} = (\mathbf{u}, \mathbf{v}, \mathbf{w})^{\top}$ is the velocity field vector, $\mathbf{p}$ is the pressure, $\boldsymbol{\tau}$ and $\boldsymbol{\tau}_t$ are respectively the viscous and turbulent stress tensors. The aim of the inverse problem here is to reconstruct the mean flow variables that are inaccurately predicted due to limitations inherent in the Boussinesq approximation

$$ \bar{\rho} \mathbf{v}^{\prime \prime} \otimes \mathbf{v}^{\prime \prime} \approx \boldsymbol{\tau}_t = \mu_t \left( \nabla \mathbf{v} + \nabla \mathbf{v}^{\top} - \frac{2}{3} (\nabla \cdot \mathbf{v}) \mathbf{I} \right) $$

used in turbulence modeling to account for the true Reynolds stress tensor $\bar{\rho} \mathbf{v}^{\prime \prime} \otimes \mathbf{v}^{\prime \prime}$. In particular, the correction inferred into the present DA-RANS formulation is introduced to compensate the modeling error in the Reynolds force vector (RFV) approximation:

$$ \nabla \cdot (\bar{\rho} \mathbf{v}^{\prime \prime} \otimes \mathbf{v}^{\prime \prime}) = \rho \mathbf{f}_{xyz} + \nabla \cdot \boldsymbol{\tau}_t, \quad (17) $$

The variables targeted by the cost function (8) usually correspond to only a subset of the conservative variables and thus do not fully describe the mean flow state, which is governed by the DA-RANS equations. Therefore, the solution of the inverse problem can reconstruct partially the exact mean flow state. Quantities not included in the measurement operator may differ from the true mean flow state. Since the DA-RANS equations must still be satisfied, the inferred correction compensates for errors in the mean flow variables that are not captured in the measurements. Consequently, the solution of the inverse problem may not be unique.

For instance, the DA-RANS momentum conservation equations simplify at non-slip solid wall boundaries as:

$$ 0 = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{f}_{xyz} \quad (18) $$

where the Reynolds stress term $\nabla \cdot \boldsymbol{\tau}_t$ and the left-hand side vanish due to the no-slip boundary condition ($\mathbf{v} = \mathbf{v}^{\prime \prime} = \mathbf{0}$). We now consider the case where the full mean velocity field is available in the cost-functional, but not the pressure field, nor the shear-stress at the wall. In such a case, the viscous strain tensor $\boldsymbol{\tau}$ should be accurately recovered by the DA procedure and an indetermination appears in the reconstruction of the wall pressure gradient: the actual pressure gradient $-\nabla p$ can be split within the sum $-\nabla p + \rho \mathbf{f}_{xyz}$, if $\mathbf{f}_{xyz}$ is left unconstrained on the wall.
Hence, to achieve a physical solution for the thermodynamic state $p$, it is therefore important to add a physical vanishing constraint on $f_{xyz}$ at the wall. At a minimum, this requires enforcing $f_{xyz} = 0$ at the wall. Preferably, regularity conditions should be applied, ensuring $f_{xyz}$ vanishes smoothly as the wall approaches.

In the following, we will enforce such conditions by constraining the values of $f_{xyz}$ to be a linear combination of gradients that exhibit such physical constraints. Rather than imposing ad hoc constraints on the tuning parameter-which are often highly dependent on the specific optimization problem-we incorporate as much information as possible from the gradient vector $g$ into the update rule of the search algorithm.

As discussed in the first part of this section, the gradient $g$ encapsulates the sensitivity of the tuning parameters with respect to the objective error, through the adjoint state. This study employs an adjoint-consistent DG scheme [23]. Specifically, along solid wall boundaries, the exact adjoint solution associated with momentum vanishes. In the above example, the adjoint consistency of the DG scheme ensures that the term $\frac{\partial R}{\partial F_{DA,k}^2}^\top \cdot Z$ in (15) is consistent with the boundary condition. Additionally, the term $\frac{\partial J_h}{\partial F_{DA,k}^2}$ also vanishes due to the norm-like properties of the penalization term, under the assumption that a homogeneous initial guess ($F_{DA}^{(0)} = 0$) is used in the optimization process. This approach allows us to maintain a general strategy that remains applicable across a wide range of tuning parameters and measurement formulations. It is therefore essential to ensure that wall forces vanish for the corrected model to be of practical relevance in adjoint-based CFD applications, such as aerodynamic shape optimization and stability-sensitivity-receptivity analyses.

## 3.2. Analysis of deterministic and stochastic optimization algorithms

The function $\mathcal{F}$ depends on the choice of optimizer. This section evaluates the limitations and advantages of common optimization algorithms. Large-scale problems, such as field-inversion in fluid mechanics, have been addressed with the L-BFGS algorithm, though the optimizer's impact has not been extensively analyzed. We begin by reviewing L-BFGS, a quasi-Newton method that is particularly useful when the Jacobian or Hessian is unavailable or expensive to compute. We then consider alternatives based on stochastic gradient descent methods, including ADAM [16] and N-ADAM, which incorporates Nesterov momentum [24]. While the suitability of N-ADAM for physics-based optimization is evaluated here, the discussion can be extended to a wide range of variations.

### 3.2.1. L-BFGS algorithm

The L-BFGS algorithm is detailed in [14,15,25]. The general update step is given by:

$$ \delta \mathbf{F}_{DA}^{(t)} = -\alpha_{W}^{(t)} \mathbf{H}^{(t)} \mathbf{g}^{(t)} \quad (19) $$

where $\alpha_{W}^{(t)}$ is the step size determined by satisfying the Wolfe conditions [26,27], and $\mathbf{H}^{(t)}$ is an approximation of the inverse Hessian matrix, using the Sherman-Morrison formula [28]. In the L-BFGS algorithm, the inverse Hessian approximation is updated as:

$$ \mathbf{H}^{(t+1)} = (\mathbf{I} - \rho_t \mathbf{s}_t \mathbf{y}_t^\top) \mathbf{H}^{(t)} (\mathbf{I} - \rho_t \mathbf{s}_t \mathbf{y}_t^\top) + \rho_t \mathbf{s}_t \mathbf{s}_t^\top $$

where $\mathbf{s}_t = \mathbf{F}_{DA}^{(t+1)} - \mathbf{F}_{DA}^{(t)}$, $\mathbf{y}_t = \nabla_{\mathbf{F}_{DA}^{(t+1)}} J_h - \nabla_{\mathbf{F}_{DA}^{(t)}} J_h$, and $\rho_t = \frac{1}{\mathbf{y}_t^\top \mathbf{s}_t}$ are used to store the last $m$ updates. At $t = 0$, the inverse Hessian matrix is initialized as a symmetric matrix:

$$ \mathbf{H}^{(0)} = \frac{\mathbf{y}_t^\top \mathbf{s}_t}{\mathbf{y}_t^\top \mathbf{y}_t} \mathbf{I}. $$

Note that, if the objective function is also convex, the initialized matrix $\mathbf{H}^{(0)}$ is also positive-definite. Thus, it is clear from (19) that the update relation $\mathcal{F}$ is linear with respect to the gradient $\mathbf{g}^{(t)}$ at each iteration.

### 3.2.2. N-ADAM algorithm

Here, we consider the N-ADAM optimizer (see Algorithms 1 and 2). Note that the analysis that follows can easily be extended to other versions of the ADAM optimizer. The corrective DOFs update (assuming $\epsilon = 0$ for simplicity), read

$$ \hat{\mathbf{m}}^{(t)} = \frac{\beta_1 \mathbf{m}^{(t-1)} + 2(1 - \beta_1) \mathbf{g}^{(t)}}{1 - \beta_1^t} \sim \mathbf{g}^{(t)} \quad (20) $$

$$ \hat{\mathbf{v}}^{(t)} = \frac{\beta_2 \mathbf{v}^{(t-1)} + 2(1 - \beta_2) \mathbf{g}^{(t)} \odot \mathbf{g}^{(t)}}{1 - \beta_2^t} \sim \mathbf{g}^{(t)} \odot \mathbf{g}^{(t)} \quad (21) $$

$$ \delta \mathbf{F}_{DA}^{(t)} = -l_r \left( \frac{\hat{\mathbf{m}}^{(t)}}{\sqrt{\hat{\mathbf{v}}^{(t)} + \epsilon}} \right) \approx -l_r \frac{\mathbf{g}^{(t)}}{\sqrt{\mathbf{g}^{(t)} \odot \mathbf{g}^{(t)}}} = -l_r \frac{\mathbf{g}^{(t)}}{|\mathbf{g}^{(t)}|} = -l_r \mathrm{sign}(\mathbf{g}^{(t)}) \quad (22) $$

where $\hat{\mathbf{m}}^{(t)}$ and $\hat{\mathbf{v}}^{(t)}$ are the first and second moments of the gradient $\mathbf{g}^{(t)}$, respectively, and $\odot$ denotes the element-wise or Hadamard product. In contrast to the L-BFGS method, the iterative update relation $\mathcal{F}$ from (22) depends non-linearly on the gradient $\mathbf{g}^{(t)}$. The non-smooth nature of its update function, combined with its sensitivity to small gradients and lack of second-order information,
can impede its performance in the specific context of field-inversion problems. The update goes to zero for very small values of the gradient, near the zero machine and assumes $\pm l_r$ values rapidly when the gradient assumes non-zero values [16]. One possible approach to mitigate this issue is the use of projected gradient descent, which constrains the solution within a predefined feasible set. However, this method often requires detailed prior knowledge of the constraints inherent to the problem and a significant effort to express them in a formal mathematical form. For this reason, in the following section, we propose a novel update rule designed to keep the solution within the physically admissible space, even when such constraints are not explicitly enforced in the optimization problem.

**Algorithm 1 N-ADAM algorithm.**

1: **given:** learning rate $l_r \in \mathbb{R}^+$, Adams exponential decay rates for moment estimates $\beta_1, \beta_2 \in [0, 1)$, Adams numerical fudge factor $\epsilon \in \mathbb{R}_{>0}$, user-defined convergence tolerances $(\epsilon_{\text{opt-tol}}, \epsilon_{\text{f-tol}}, \epsilon_{\text{g-tol}}) \in \mathbb{R}^3_{>0}$

2: **given:** Objective function $J_h$

3: **initialize:** time step $t \leftarrow 0$, solution $U^{(t=0)}$ and parameter $F_{\text{DA}}^{(t=0)}$ DOFs, first moment vector $\mathbf{m}^{(t=0)} \leftarrow 0$, second moment vector $\mathbf{v}^{(t=0)} \leftarrow \mathbf{0}$

4: **while** $F_{\text{DA}}^{(t)}$ is not converged **do**

5: $t \leftarrow t + 1$

6: $\nabla_{F_{\text{DA}}^{(t)}} J_h \leftarrow \text{AdjointSolver}(U^{(t-1)}, F_{\text{DA}}^{(t-1)})$

▷ Eqs. (11) and (15)

7: $\mathbf{g}^{(t)} \leftarrow \nabla_{F_{\text{DA}}^{(t)}} J_h$

8: $\mathbf{m}^{(t)} \leftarrow \beta_1 \mathbf{m}^{(t-1)} + (1 - \beta_1) \mathbf{g}^{(t)}$

9: $\mathbf{v}^{(t)} \leftarrow \beta_2 \mathbf{v}^{(t-1)} + (1 - \beta_2) \mathbf{g}^{(t)} \odot \mathbf{g}^{(t)}$

10: $\hat{\mathbf{m}}^{(t)} \leftarrow \frac{\mathbf{m}^{(t)} + (1 - \beta_1) \mathbf{g}^{(t)}}{1 - \beta_1^t}$

11: $\hat{\mathbf{v}}^{(t)} \leftarrow \frac{\mathbf{v}^{(t)}}{1 - \beta_2^t}$

12: $\mathbf{F}_{\text{DA}}^{(t)} = \mathbf{F}_{\text{DA}}^{(t-1)} - l_r \left( \frac{\hat{\mathbf{m}}^{(t)}}{\sqrt{\hat{\mathbf{v}}^{(t)} + \epsilon}} \right)$

▷ Update DOFs

13: $U^{(t)}, J_h^{(t)} \leftarrow \text{DirectSolver}(U^{(t-1)}, F_{\text{DA}}^{(t)})$

▷ Eqs. (10) and (8)

14: **if** CHECKCONVERGENCE($J_h^{(t)}, J_h^{(t-1)}, \nabla_{F_{\text{DA}}^{(t)}} J_h, \epsilon_{\text{opt-tol}}, \epsilon_{\text{f-tol}}, \epsilon_{\text{g-tol}}$) **then**

15: **break**

16: **end if**

17: **end while**

18: **return:** Optimized parameters of DOFs $F_{\text{DA}}^{(t)}$

**Algorithm 2 Procedure to evaluate convergence criteria of $J_h$ minimization.**

1: **function** CHECKCONVERGENCE($J_h^{(t)}, J_h^{(t-1)}, \nabla_{F_{\text{DA}}^{(t)}} J_h, \epsilon_{\text{opt-tol}}, \epsilon_{\text{f-tol}}, \epsilon_{\text{g-tol}}$)

2: **if** $J_h^{(t)} \le \epsilon_{\text{opt-tol}}$ **then**

▷ First criterion: absolute value of objective

3: **return true**

4: **else if** $|J_h^{(t)} - J_h^{(t-1)}| \le \epsilon_{\text{f-tol}} \max(1.0, J_h^{(t)}, J_h^{(t-1)})$ **then**

▷ Second criterion: relative change in objective

5: **return true**

6: **else if** $\|\nabla_{F_{\text{DA}}^{(t)}} J_h\|_{\infty} \le \epsilon_{\text{g-tol}}$ **then**

▷ Third criterion: gradient norm

7: **return true**

8: **end if**

9: **return false**

10: **end function**

### 3.2.3. V-N-ADAM-DG algorithm

This section introduces a new stochastic optimizer, termed V-N-ADAM-DG, which incorporates gradient knowledge in a similar way such as L-BFGS. Similar version has been proposed by [17] for geometric optimization problems where vectorial quantities are optimized. Here, a further adaptation to the DG representation of the corrective term is taken into account. Its characteristics is that each component and each $l$th DOF of the corrective DOFs vector $(F_{\text{DA},\kappa}^l)_{1 \le l \le N_p} = (F_{1,\kappa}^l, \dots, F_{N_r,\kappa}^l)_{1 \le l \le N_p}$ is updated separately into the optimization procedure (see Algorithm 3).

The gradient of the objective function with respect to the corrective DOFs, $\nabla_{F_{\text{DA}}^{(t)}} J_h$, initially has dimensions $N_{\text{DOFs}} \times N_r$, where $N_r$ denotes the number of scalar components in the corrective DOFs vector, and $N_{\text{DOFs}} = N_{\text{elts}} \times N_p$ represents the total number of DOFs per component. Each column vector $\mathbf{g}_r$ of the gradient is used to update the respective component of the corrective DOFs vector
**Algorithm 3 V-N-ADAM-DG algorithm.**

1: **given:** learning rate $l_r \in \mathbb{R}^+$, Adams exponential decay rates for moment estimates $\beta_1, \beta_2 \in [0, 1)$, Adams numerical fudge factor $\epsilon \in \mathbb{R}_{>0}$, user-defined convergence tolerances $(\epsilon_{\text{opt-tol}}, \epsilon_{\text{f-tol}}, \epsilon_{\text{g-tol}}) \in \mathbb{R}^3$

2: **given:** Number of components $N_r$ of $\mathbf{F}_{\text{DA}}$, number of total DOFs $N_{\text{DOFs}} = N_{\text{elts}} \times N_p$

3: **given:** Objective function $J_h$

4: **initialize:** time step $t \leftarrow 0$, solution $\mathbf{U}^{(t=0)}$ and parameter $\mathbf{F}_{\text{DA}}^{(t=0)}$ DOFs

5: **initialize:** first moment matrix $\mathbf{m}^{(t=0)} \leftarrow 0$ of shape $N_{\text{elts}} \times N_r^{\text{DG}}$ with $N_r^{\text{DG}} = N_r \times N_p$

6: **initialize:** second moment vector $\mathbf{v}^{(t=0)} \leftarrow 0$ of shape $N_r^{\text{DG}} \times 1$

7: **while** $\mathbf{F}_{\text{DA}}^{(t)}$ is not converged **do**

8: &nbsp;&nbsp;&nbsp;&nbsp; $t \leftarrow t + 1$

9: &nbsp;&nbsp;&nbsp;&nbsp; $\nabla_{\mathbf{F}_{\text{DA}}^{(t)}} J_h \leftarrow \text{AdjointSolver}(\mathbf{U}^{(t-1)}, \mathbf{F}_{\text{DA}}^{(t-1)})$

▷ Eqs. (10)-(11) and (15). Gradient has shape $N_{\text{DOFs}} \times N_r$

10: **for** $r \in [1, N_r]$ **do**

11: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{g}_r^{(t)} \leftarrow \nabla_{\mathbf{F}_{\text{DA},r}} J_m$

▷ Gradient of $J_h$ w.r.t. the $r$th component of $\mathbf{F}_{\text{DA}}^{(t)}$

12: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{G}_r^{(t)} \leftarrow \mathbf{g}_r^{(t)}$

▷ Reshape gradients from $N_{\text{DOFs}} \times 1$ to $N_{\text{elt}} \times N_p$

13: **for** $l \in [1, N_p]$ **do**

14: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{G}_{r,l}^{(t)} \leftarrow \mathbf{G}_r^{(t)} = (\mathbf{G}_{r,l}^{(t)})_{1 \le l \le N_p}$

▷ Gradient w.r.t the $l$th DOF of $\mathbf{F}_{\text{DA},r}$

15: &nbsp;&nbsp;&nbsp;&nbsp; $k \leftarrow (r-1)N_p + l$

16: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{m}_k^{(t)} \leftarrow \beta_1 \mathbf{m}_k^{(t-1)} + (1-\beta_1) \mathbf{G}_{r,l}^{(t)}$

17: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{v}_k^{(t)} \leftarrow \beta_2 \mathbf{v}_k^{(t-1)} + (1-\beta_2) \|\mathbf{G}_{r,l}^{(t)}\|_2^2$

18: &nbsp;&nbsp;&nbsp;&nbsp; $\hat{\mathbf{m}}_k^{(t)} \leftarrow \frac{\mathbf{m}_k^{(t)} + (1-\beta_1)\mathbf{G}_{r,l}^{(t)}}{1-\beta_1^t}$

19: &nbsp;&nbsp;&nbsp;&nbsp; $\hat{v}_k^{(t)} \leftarrow \frac{v_k^{(t)}}{1-\beta_2^t}$

$$20: \quad \mathbf{F}_{\text{DA},r}^{l,(t)} = \mathbf{F}_{\text{DA},r}^{l,(t-1)} - l_r \left( \frac{\hat{\mathbf{m}}_k^{(t)}}{\sqrt{\hat{v}_k^{(t)} + \epsilon}} \right)$$

▷ Update DOFs

21: **end for**

22: **end for**

23: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{U}^{(t)}, J_h^{(t)} \leftarrow \text{DirectSolver}(\mathbf{U}^{(t-1)}, \mathbf{F}_{\text{DA}}^{(t)})$

▷ Eqs. (10) and (8)

24: **if** CHECKCONVERGENCE($J_h^{(t)}, J_h^{(t-1)}, \nabla_{\mathbf{F}_{\text{DA}}^{(t)}} J_h, \epsilon_{\text{opt-tol}}, \epsilon_{\text{f-tol}}, \epsilon_{\text{g-tol}}$) **then**

25: **break**

▷ Convergence criteria from Algorithm 2

26: **end if**

27: **end while**

28: **return:** Optimized parameters of DOFs $\mathbf{F}_{\text{DA}}^{(t)}$

independently (as shown in line 10 of Algorithm 3). Specifically, $\mathbf{g}_r$ denotes the gradient of the objective function with respect to the $r$th component of the corrective parameter $(\mathbf{F}_{\text{DA},r}^{l})_{1 \le l \le N_p}$.

This gradient is reshaped into a $N_{\text{elts}} \times N_p$ matrix $\mathbf{G}_r$, and each column is updated through a second loop over $N_p$. This structure preserves the scaled modal coefficients used in polynomial reconstruction (see line 13 of Algorithm 3). Thus, the $l$th DOF of $(\mathbf{F}_{\text{DA},r}^{l})_{1 \le l \le N_p}$ is updated using $\mathbf{G}_{r,l}$ exclusively, which corresponds to the gradient of the $l$th DOF.

Once the original gradient $\nabla_{\mathbf{F}_{\text{DA}}^{l}} J_h$ is reshaped, the final modification lies in the definition of the second momentum: the element-wise product of the gradient used in traditional ADAM-based algorithms (line 9 of Algorithm 1) is replaced by the Euclidean norm (line 17 of Algorithm 3). Consequently, the sizes of the first and second momentum terms become $N_{\text{elts}} \times N_r^{\text{DG}}$ and $N_r^{\text{DG}} \times 1$, respectively, with $N_r^{\text{DG}} = N_r \times N_p$.

For each $r$th component of $\mathbf{F}_{\text{DA}}$ and each $l$th DOF, the initial update reads:

$$\hat{\mathbf{m}}_k^{(t)} = \frac{\beta_1 \mathbf{m}_k^{(t-1)} + 2(1-\beta_1)\mathbf{G}_{r,l}^{(t)}}{1-\beta_1^t} \sim \mathbf{G}_{r,l}^{(t)} \quad (23)$$

$$\hat{v}_k^{(t)} = \frac{\beta_2 \mathbf{v}_k^{(t-1)} + 2(1-\beta_2) \|\mathbf{G}_{r,l}^{(t)}\|_2^2}{1-\beta_2^t} \sim \|\mathbf{G}_{r,l}^{(t)}\|_2^2 \quad (24)$$

$$\delta \mathbf{F}_{\text{DA},r}^{l,(t)} = -l_r \left( \frac{\hat{\mathbf{m}}_k^{(t)}}{\sqrt{\hat{v}_k^{(t)} + \epsilon}} \right) \approx -l_r \frac{\mathbf{G}_{r,l}^{(t)}}{\|\mathbf{G}_{r,l}^{(t)}\|_2} \quad (25)$$
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Update rule visualization at a general *l* iteration: Comparison of the update rules for the N-ADAM (blue arrow) and V-N-ADAM-DG (red arrow) optimizers in the case of $N_p = 1$, considering the $N_r = 2$ components of the optimal parameter $\mathbf{F}_{\text{DA}} = (F_x, F_y)$. The inverse gradient direction (descent direction) is represented by the black arrow.</figcaption></figure>

Here, the update function $\mathbf{F}$ is linear, as a result of the modifications introduced. The update for each rth component and each lth DOF is linear with respect to the normalized gradient $\mathbf{G}_{\text{r},l}^{(t)}$ at each iteration. As a result, the update converges linearly to zero as the gradient diminishes. The optimization process is also regularized by the normalization introduced through the second momentum, which constraints the update in a bounded $[-l_r, +l_r]$ range. A graphical representation of the key differences between the two stochastic optimizers is presented in Fig. 1 for the case of $N_p = 1$, considering only the $N_r = 2$ components of the optimal parameter (referred as $\mathbf{F}_{\text{DA}} = (F_x, F_y)$ in the scheme). The update direction of the N-ADAM optimizer (22) is misaligned with the true gradient direction, causing the parameter updates to follow diagonal paths. In contrast, the modified update rule (25) ensures that updates follow the gradient direction while using the gradient norm in the measurement space to adapt the update intensity. This approach prevents excessive or intrusive updates in regions with weak gradients, improving overall stability and accuracy.

### 3.3. Invariance under linear transformations and rotation of the physical reference system

A suitable property of an optimizer is its invariance under non-singular linear transformations, i.e. transformations between vector spaces that are both linear and invertible. As detailed in [29], quasi-Newton methods, such as L-BFGS, are known to exhibit this property. Consider the nth component of the corrective DOFs $\mathbf{F}_{\text{DA}}$ (with $1 \le n \le N_{\text{DOFs}}$) being optimized. When a non-singular transformation is applied:

$$ \hat{\mathbf{F}}_{\text{DA}} = \mathbf{A}\mathbf{F}_{\text{DA}} + \mathbf{B} \qquad (26) $$

where $\mathbf{A} \in \mathbb{R}^{N_r \times N_r}$ is invertible and $\mathbf{B} \in \mathbb{R}^{N_r}$.

From a physical standpoint, the corrective terms to be optimized can be defined as vectorial quantities and this transformation could correspond to the specific case of a reference system rotation, represented by

$$ \mathbf{Q} = \begin{bmatrix} \cos \alpha_y \cos \alpha_z & \sin \alpha_x \sin \alpha_y \cos \alpha_z - \cos \alpha_x \sin \alpha_z & \cos \alpha_x \sin \alpha_y \cos \alpha_z + \sin \alpha_x \sin \alpha_z \\ \cos \alpha_x \sin \alpha_z & \sin \alpha_x \sin \alpha_y \sin \alpha_z + \cos \alpha_x \cos \alpha_z & \cos \alpha_x \sin \alpha_y \sin \alpha_z - \sin \alpha_x \cos \alpha_z \\ -\sin \alpha_y & \sin \alpha_x \cos \alpha_y & \cos \alpha_x \cos \alpha_y \end{bmatrix} $$

where any rotation around the x, y and z axes is quantified by an angle $\alpha_x, \alpha_y, \alpha_z$, respectively. In the DA-RANS problem, the corrective term $\mathbf{f}_{\text{DA}}$ incorporates mixed rotational properties since, in the baseline RANS equations, the momentum equations rotate with the reference system, whereas the energy equation is rotation-invariant (the properties are inherited from the general Navier-Stokes equations, where additional proofs are referred to [30]). Given a rotation $\mathbf{Q}$, such that $\mathbf{Q}^\mathrm{T}\mathbf{Q} = \mathbf{I}$, in the momentum equation:

$$ \rho \mathbf{f}_{xyz}^{\mathbf{Q}} + (\nabla \cdot \boldsymbol{\tau}_i)^{\mathbf{Q}} = \rho \mathbf{f}_{xyz}^{\mathbf{Q}} + \mathbf{Q} \nabla \cdot \boldsymbol{\tau}_i, \qquad (27) $$

the corrective term must rotate consistently with the rest of the equations, i.e.,

$$ \mathbf{f}_{xyz}^{\mathbf{Q}} = \mathbf{Q}\mathbf{f}_{xyz}. \qquad (28) $$

to avoid altering the momentum equation.

By storing the DOFs of the inferred correction in the vector $\mathbf{F}_{\text{DA}} = (\mathbf{F}_{\text{DA}}^l)_{1 \le l \le N_p}$, we can define a new rotation matrix $\tilde{\mathbf{Q}}$, which applies the standard rotation over the three axis for each modal DOF $l$ as:

$$ (\mathbf{F}_{\text{DA}})^{\tilde{\mathbf{Q}}} = (\mathbf{Q}\mathbf{F}_{\text{DA}}^{l})_{1 \le l \le N_p} = \tilde{\mathbf{Q}}\mathbf{F}_{\text{DA}}, \qquad (29) $$

which is a specific case of transformation (26), where $\mathbf{A} = \tilde{\mathbf{Q}}, \mathbf{B} = \mathbf{0}$ and $N_r = 3$.

Finally, under the assumption that the objective functions remain invariant under the given transformation, it can be demonstrated that the gradient $\mathbf{g}$ exhibits equivariance with respect to the new transformation matrix, such that $\mathbf{g}^{\tilde{\mathbf{Q}}} = \tilde{\mathbf{Q}}\mathbf{g}$. However, the N-ADAM
update rule (22) undergoes the following transformation:

$$ (\delta \mathbf{F}_{\text{DA}}^{(t)})^{\tilde{Q}} = -l_r \operatorname{sign}(\tilde{\mathbf{Q}}\mathbf{g}^{(t)}) \neq \tilde{\mathbf{Q}}\delta\mathbf{F}_{\text{DA}}^{(t)}, \qquad (30) $$

whereas the V-N-ADAM-DG update rule (25) transforms as:

$$ (\delta \mathbf{F}_{\text{DA},r}^{l,(t)})^{\tilde{Q}} = -l_r \frac{\tilde{\mathbf{Q}}\mathbf{G}_{\text{r},l}^{(t)}}{\|\tilde{\mathbf{Q}}\mathbf{G}_{\text{r},l}^{(t)}\|_2^2} = -l_r \frac{\tilde{\mathbf{Q}}\mathbf{G}_{\text{r},l}^{(t)}}{\|\mathbf{G}_{\text{r},l}^{(t)}\|_2^2} = \tilde{\mathbf{Q}}\delta\mathbf{F}_{\text{DA},r}^{l,(t)}. \qquad (31) $$

While full invariance under the general transformation in (26) is not achieved, the V-N-ADAM update rule's invariance under transformations of the form $f_{d-\text{rot}}: \mathbf{x} \to \mathbf{Q}\mathbf{x}$ (or equivalently $f_{(d+1)-\text{rot}}: \mathbf{F}_{\text{DA}} \to \tilde{\mathbf{Q}}\mathbf{F}_{\text{DA}}$) offers a significant advantage. It ensures that the optimizer's search for the optimal solution is unaffected by the choice of reference system, making it more robust.

### 3.4. Computational time and memory requirements estimation

The three algorithms are evaluated here in terms of computational time. The total time for the DA procedure, $T_{\text{DA}}$, can be decomposed into three main contributions $T_{\text{DA}} = T_{\text{U}} + T_{\text{Z}} + T_{\text{opt}}$, where $T_{\text{U}}$ represents the time required to solve the direct problem (10) and evaluate the objective function via eq. (8), $T_{\text{Z}}$ corresponds to the time for solving the adjoint problem (11) and computing the gradient through (15), and $T_{\text{opt}}$ accounts for the time to update the DOFs according to the search algorithm.

The first two contributions depend primarily on the capabilities of the numerical solver coupled with the optimizer, while $T_{\text{opt}}$ depends solely on the definition of the update rule. As detailed in [14], the time per iteration of the L-BFGS algorithm is approximately $(4m+1)N_{\text{DOFs}}^{\text{DA}}$, where the number of stored updates for the Hessian approximation is typically $3 \le m \le 20$. No explicit estimates of the iteration time for the N-ADAM algorithm are available in the literature; however, $T_{\text{opt}}$ for both N-ADAM and V-N-ADAM-DG can be inferred from Algorithms 1 and 3 and also scales as $\mathcal{O}(N_{\text{DOFs}}^{\text{DA}})$. In practice, we observed a cost $cN_{\text{DOFs}}^{\text{DA}}$ with $10 \le c \le 20$ depending on implementation details. Hence, the computational times of these optimizers are likely comparable.

While $T_{\text{U}}$ and $T_{\text{Z}}$ can be marginally influenced by the choice of search algorithm, their contribution is generally dominated by the numerical parameters set in the solver (e.g., the strategy used to solve linear systems). This discussion naturally extends to memory requirements. For the direct and adjoint problems, memory demands are fixed since they are treated independently of the search algorithm. L-BFGS requires roughly $\mathcal{O}(mN_{\text{DOFs}}^{\text{DA}})$ storage locations, which is already efficient for very large-scale problems. In contrast, N-ADAM and V-N-ADAM-DG only need to store a single previous update for the Nesterov momentum acceleration, resulting in memory requirements of order $\mathcal{O}(N_{\text{DOFs}}^{\text{DA}})$, slightly lower than those of L-BFGS.

The overall computational cost of the optimization problem depends on the total number of iterations required for convergence. In this respect, L-BFGS is advantageous, as the use of curvature information typically leads to fewer iterations. In the present work, V-N-ADAM-DG has been designed without specific attention to convergence acceleration. Consequently, it inherits the convergence speed of N-ADAM, while differing from it in the quality of the obtained optimal solution. Although V-N-ADAM-DG may be slower than L-BFGS, it is capable of handling stiff cost functions that involve a combination of different metrics, making it more robust in such scenarios. Combination with adaptive learning rate strategies and a detailed convergence study could be the subject of future work.

## 4. Application to mean flow reconstruction in RANS simulations discretized using a high-order method

### 4.1. Baseline RANS equations

In this study, we evaluate the findings of the present work in the context of reconstructing mean flow variables in RANS simulations. The focus is on the resolution of field-inversion problems for the compressible RANS equations, discretized using a DG formalism. The baseline RANS model is the Spalart-Allmaras (SA) one-equation turbulence model [31], with the correction proposed in [32] to account for shear flows at low Reynolds numbers.

Let $\Omega \in \mathbb{R}^d$ be a bounded domain, where $d$ denotes the spatial dimension. The baseline steady-state SA-RANS equations for the conservative variables $\mathbf{u} = (\rho, \rho\mathbf{v}^T, \rho\mathbf{E}, \rho\tilde{\mathbf{v}})^T$ are expressed in vector form as:

$$ \nabla \cdot \mathbf{f}_c(\mathbf{u}) - \nabla \cdot \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) \qquad (32) $$

in $\Omega$, subject to prescribed boundary conditions on $\partial\Omega$.

In eq. (32), the nonlinear convective flux $\mathbf{f}_c(\mathbf{u})$, the viscous flux $\mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u})$, and the source term $\mathbf{s}(\mathbf{u}, \nabla \mathbf{u})$ are given by:

$$ \mathbf{f}_c(\mathbf{u}) = \begin{pmatrix} \rho \mathbf{v}^T \\ \rho \mathbf{v} \mathbf{v}^T + \mathbf{p} \mathbf{I} \\ \rho \mathbf{H} \mathbf{v}^T \\ \rho \tilde{\mathbf{v}} \mathbf{v}^T \end{pmatrix}, \quad \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \begin{pmatrix} 0 \\ \tau + \tau^T \\ \mathbf{v}^T (\tau + \tau^T) - \mathbf{q}^T - \mathbf{q}^T \\ \frac{1}{\sigma} \zeta \nabla \tilde{\mathbf{v}}^T \end{pmatrix}, \quad \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) = \begin{pmatrix} 0 \\ 0 \\ 0 \\ c_{b1} \hat{S} \rho \tilde{\mathbf{v}} - c_{w1} f_w \frac{\rho \tilde{\mathbf{v}}^2}{d_w^2} + \frac{c_{b2} \rho}{\sigma} \nabla \tilde{\mathbf{v}} \cdot \nabla \tilde{\mathbf{v}} \end{pmatrix}. $$

All terms of the SA model are provided in Appendix A.
## 4.2. DA-RANS equations: Reynolds force correction

This section introduces the first formulation considered for the field-inversion problem. Previous studies have successfully applied a similar methodology to incompressible flows [20]. The incorporation of source terms in the momentum equations has demonstrated strong performance across various configurations, including turbulent flow over a square cylinder and flow over periodically arranged hills in a channel [11,33]. In [13], this methodology was extended to the compressible RANS Eq. (32) by introducing a vector source term in the momentum and energy equations:

$$ \nabla \cdot \mathbf{f}_c(\mathbf{u}) - \nabla \cdot \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) + \rho \mathbf{f}_{DA}, \quad (33) $$

where

$$ \mathbf{f}_{DA} = \begin{pmatrix} 0 \\ \mathbf{f}_{xyz} \\ f_e + \mathbf{v}^\top \mathbf{f}_{xyz} \\ 0 \end{pmatrix}. $$

The correction term $\mathbf{f}_{DA}$ is sought within the DG function space, i.e., $\mathbf{f}_{DA} \in (\mathcal{V}_h^p)^{d+3}$, with $\mathbf{f}_{xyz} = (f_x, f_y, f_z)^\top$ and:

$$ f_x = \sum_{l=1}^{N_p} \phi_k^l(\mathbf{x}) F_{x,\kappa}^l, \quad f_y = \sum_{l=1}^{N_p} \phi_k^l(\mathbf{x}) F_{y,\kappa}^l, \quad f_z = \sum_{l=1}^{N_p} \phi_k^l(\mathbf{x}) F_{z,\kappa}^l, \quad f_e = \sum_{l=1}^{N_p} \phi_k^l(\mathbf{x}) F_{e,\kappa}^l. $$

The DOFs associated with the representation of $\mathbf{f}_{DA}$ within the element $\kappa$ are organized in the vector

$$ \mathbf{F}_{DA,\kappa} = (F_{x,\kappa}^l, F_{y,\kappa}^l, F_{z,\kappa}^l, F_{e,\kappa}^l)_{1 \le l \le N_p}^\top. $$

The total number of DOFs for $\mathbf{F}_{DA}$ is given by $N_{DOFs} = N_{elts} \times N_p \times N_r$, where $N_{elts}$ denotes the number of elements in $\Omega_h$, and $N_r = d + 1$ represents the number of components of the corrective parameter.

## 4.3. DA-RANS equations: Turbulent transport equation correction

The second formulation introduces an additional source term $\rho \tilde{v} f_{\tilde{v}}$ in the SA equation, as proposed by Cato et al. [11]:

$$ \nabla \cdot \mathbf{f}_c(\mathbf{u}) - \nabla \cdot \mathbf{f}_v(\mathbf{u}, \nabla \mathbf{u}) = \mathbf{s}(\mathbf{u}, \nabla \mathbf{u}) + \rho \mathbf{f}_{DA}, \quad (34) $$

where

$$ \mathbf{f}_{DA} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ \tilde{v} f_{\tilde{v}} \end{pmatrix}. $$

The inferred correction function $f_{\tilde{v}}$ is designed to modify the eddy-viscosity values by adjusting the source terms of the SA turbulence model equation. The correction term $f_{\tilde{v}}$ is pre-multiplied by $\rho \tilde{v}$ to spatially localize its effect within regions where $\rho \tilde{v}$ is large, corresponding to areas of significant turbulence activity. This source term is sought within the DG function space $\mathcal{V}_h^p$:

$$ f_{\tilde{v}}(\mathbf{x}) = \sum_{l=1}^{N_p} \phi_k^l(\mathbf{x}) F_{\tilde{v},\kappa}^l \quad \forall \mathbf{x} \in \kappa, \quad \kappa \in \Omega_h, \quad (35) $$

where $\mathbf{F}_{DA,\kappa} = (F_{\tilde{v},\kappa}^l)_{1 \le l \le N_p}$ represents the DOFs associated with $f_{\tilde{v}}(\mathbf{x})$ in the element $\kappa$. Given that the corrective parameter consists of a single component, i.e., $N_r = 1$, the total number of DOFs for $\mathbf{F}_{\tilde{v}}$ is given by $N_{DOFs} = N_{elts} \times N_p$.

# 5. Numerical results

## 5.1. Test case configurations

In this section, we apply the proposed field-inversion strategies to two test cases of increasing complexity. Numerical results for a laminar flow configuration (referred to as L-CYL) are used to assess the ability of the optimizers to solve inverse problems with high-order polynomial degrees in the DG framework. A supersonic shock-wave/turbulent boundary layer interaction configuration (referred to as SWTBLI) is then investigated to show the response of the optimizers to treat inverse problems with different cost functions and inferred corrections. Although both configurations are two-dimensional, they provide representative and well-controlled scenarios for assessing the proposed field-inversion methodology and the optimization algorithms. The L-CYL case emphasizes the influence of numerical discretization on the optimal solution and optimizer behavior, while the SWTBLI case represents a more complex flow regime featuring strong shear-layer separation, shock waves, and turbulence-modeling challenges, thereby paving the way for extending the present methodology to realistic aerodynamic applications.
<figure><img src="image_4.png" /><figcaption>Fig. 2. L-CYL (error on velocity field): Mean streamwise momentum contour $\frac{\rho u}{\rho_{\infty} u_{\infty}} \in [-0.1, 1.1]$. Black continuous lines refer to the DNS solution and contours refer to the baseline and DA-RANS solutions. (a)Baseline; (b)L-BFGS; (c)N-ADAM; (d)V-N-ADAM-DG.</figcaption></figure>

The numerical results are obtained with the Aghora CFD code, a high-order DG solver for compressible turbulent flow simulations developed at ONERA [34]. The RANS equations were sought in non-dimensional form and normalization was done with free-stream quantities ($\rho_{\infty}$, $v_{\infty}$, $T_{\infty}$) and a reference length scale (the diameter of the cylinder for the L-CYL case and the inlet boundary layer thickness for the SWTBLI configuration). The solutions of the baseline (2) and DA (3) problems are obtained via a quasi-Newton method, incorporating a pseudo-time stepping approach with local time stepping. Given the sparsity of the implicit system matrix and the large number of DOFs, a Krylov subspace method is employed for solving the linear system at each nonlinear iteration. Specifically, we utilize a restarted GMRES method [35] with a block ILU(0) preconditioner [36]. To enhance computational efficiency, a Jacobian-free strategy is adopted, as detailed in [13].

## 5.2. Laminar flow past a cylinder (L-CYL)

We begin by analyzing the flow around a cylinder (L-CYL) at a low Reynolds number $Re = 80$ and Mach number $M = 0.15$, based on the cylinder diameter and freestream parameters. In this scenario, the flow is unsteady and laminar, so this test serves solely for validation purposes, to verify that the optimizer can reproduce the time-averaged solution of the unsteady laminar flow. The computational grid consists of $N_{\text{elt}} = 4386$ cubic quadrangles and extends up to 30 cylinder diameters. No-slip adiabatic boundary conditions are applied at the wall, and non-reflecting boundary conditions are enforced at the artificial boundaries. In this case, $\Omega_m = \Omega_h$, and the reference data were obtained from a time-averaged DNS on the same mesh. A baseline computation is carried out using a polynomial reconstruction order $p = 2$. The contour plot of $\frac{\rho u}{\rho_{\infty} u_{\infty}}$ for the baseline solution is presented in Fig. 2a and compared to the DNS (black isolines). The high-order spatial discretization, achieved through the use of high polynomial degrees, enables a satisfactory representation of the solution on a coarse mesh. Upon comparison of the baseline and DNS results, it becomes evident that the baseline simulation overestimates the wake size.

For the field-inversion problem (7), the mean velocity field is included in the error metric:

$$J(v) = \frac{1}{2} \|v_h - v_m\|_{L^2(\Omega_h)}^2 \qquad (36)$$
<figure><img src="image_5.png" /><figcaption>Fig. 3. L-CYL (error on velocity field): Evolution of the normalized cost function for the L-BFGS (blue squares), N-ADAM (green deltas) and V-N-ADAM-DG (orange left triangles) optimizers.</figcaption></figure>

Significant correction is achieved with each optimizer, and the DA-RANS solutions are compared to DNS results in Fig. 2b to d. However, when analyzing the objective function history (see Fig. 3), only L-BFGS and V-N-ADAM-DG reach lower minima. The N-ADAM optimizer shows less stable convergence compared to the others, and the optimal corrective forcing term in the momentum equations deviates significantly from the L-BFGS and V-N-ADAM optima (Fig. 4). The temporary increase in the objective function observed during the final iterations may be attributed to the stochastic nature of the optimizer, which explores various regions within the solution space. This fluctuation suggests that the V-N-ADAM-DG algorithm examined a region with higher objective values before identifying a more favorable direction toward convergence at the L-BFGS minimum.

In Fig. 4c and d, a non-zero forcing term is observed at the solid wall boundary, indicating that the forcing term from the N-ADAM optimum does not adequately represent the RFV correction. In contrast, as shown in Fig. 4e and f, the forcing term smoothly approaches zero near the solid wall, with the correction concentrated in the wake region downstream of the cylinder.

An additional $L^2$-norm-based penalty on the inferred correction is tested in the present problem to regularize the optimal solution space. For L-BFGS, this penalty is implemented as a standard $L^2$-regularization term. In contrast, for N-ADAM and V-N-ADAM-DG, the penalty is incorporated using the decoupled regularization strategy proposed by [37], in which the gradient of the penalty term is treated separately from the gradient associated with the measurement errors. This formulation corresponds to the standard implementation of the N-ADAM algorithm in the literature and has been extended to the novel optimizer. Several levels of penalty intensity, tuned by the parameter $\alpha_p$ in Eq. (8), have been examined. The N-ADAM optimal solution remains unchanged compared to the non-regularized distribution shown in Fig. 4c and d, whereas both L-BFGS and V-N-ADAM-DG exhibit improvements. In these cases, the penalty effectively suppresses excessive correction amplitudes without significantly increasing the measurement error.

The optimal solutions obtained by L-BFGS and V-N-ADAM-DG are consistent with the results found by Foures et al. (2014) [38]. However, in their work, the open-source FreeFEM solver [39] based on finite elements methods was used for the corrective term, while here, a high-order DG method is employed. The use of high-order polynomial reconstruction increases the complexity of the optimization problem, as it expands the size of the optimization space. Additionally, the DG modal basis is hierarchical, meaning that higher-order modes can be added without altering the lower-order basis. As a result, the primary corrective degrees of freedom (DOFs) should correspond to the lower-order modes, with higher-order modes serving to enhance the accuracy of the solution representation. Consequently, it is expected that the DOF coefficients decay rapidly, as demonstrated by [36] in the DG context and by [40] for analytical functions expressed in a modal Legendre basis. To assess this, a DOFs decay factor can be defined as

$$ \lambda(F_r^l) = \frac{\int_{\Omega_h} F_{r,k}^l d\mathbf{x}}{\int_{\Omega_h} F_{r,k}^{l=1} d\mathbf{x}} $$

The decay factor for each r component of $F_{DA}$ is shown as a function of the degrees of the modal basis in Fig. 5. The L-BFGS and V-N-ADAM-DG optimizers produce better DOF decay, whereas the N-ADAM optimum gives nearly equal weight to all DOFs in the polynomial reconstruction, making the optimal form highly sensitive to the polynomial degree p. The absence of decay in the DOFs obtained with N-ADAM as l increases indicates that the optimal solution, $f_{DA}$, exhibits a loss of regularity, as previously shown in Fig. 4c and d. This lack of decay in the N-ADAM solution DOFs can be attributed to the algorithm's update rule. Specifically, the N-ADAM update rule allows DOFs to adopt values of $\pm l_r$ independently of the modal degree l they represent. By updating each of the $N_p$ DOFs at the same rate, regardless of the degree l in the polynomial representation, N-ADAM promotes an equal contribution across all DOFs in the search for the optimum. In contrast, V-N-ADAM-DG addresses this issue through an update rule that allows each DOF to vary within a continuous range $[+l_r, -l_r]$, permitting more nuanced adjustments. Furthermore, in V-N-ADAM-DG update rule, the lth DOF is adjusted based only on the gradient of $J_h$ with respect to the lth DOF of $F_{DA}$. This enables a more tailored optimization process and promotes smoother corrections, as greater sensitivity of $J_h$ is assigned to the lower-order modal coefficients.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** L-CYL (error on velocity field): Streamwise and vertical corrective forces colored contour and stream-wise velocity black isolines for the L-BFGS (a,b), N-ADAM (c,d) and V-N-ADAM-DG (e,f) optimizers. (a)ρfₓ (L-BFGS); (b)ρfᵧ (L-BFGS); (c)ρfₓ (N-ADAM); (d)ρfᵧ (N-ADAM); (e)ρfₓ (V-N-ADAM-DG); (f)ρfᵧ (V-N-ADAM-DG).</figcaption></figure>

5.3. Shock-wave/turbulent boundary layer interaction (SWTBLI): Momentum and energy correction terms

In this section, we consider the complex problem of shock-wave/turbulent boundary layer interaction (SWTBLI). When a shock
wave interacts with a boundary layer, it triggers intricate phenomena, including the deceleration of the boundary layer flow and the
propagation of the shock wave in a multi-layered structure, as detailed in [41–43]. This configuration has been widely investigated
through experimental studies [44,45] and numerical simulations [46–48].

The current analysis focuses on a case with an inlet Mach number of M∞ = 2.28, adiabatic walls, and a deflection angle φ = 9.5°,
based on the DNS database from Volpiani et al. [47]. For the RANS simulations, a computational domain Ωh = (7,55)δ0 × (0,15)δ0 is
used, discretized with a mesh of Nₓ × Nᵧ = 350 × 250 elements. A polynomial degree of p = 1 is applied for the DG representation,
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** L-CYL (error on velocity field): DOFs decay factor for the streamwise momentum (a), vertical momentum (b) and energy (c) components. Comparison of the L-BFGS (blue squares), N-ADAM (green deltas) and V-N-ADAM-DG (orange left triangles) optimizers for $1 \le l \le N_p = 6$ with $p = 2$. (a) $\lambda(F_x^l)$; (b) $\lambda(F_y^l)$; (c) $\lambda(F_e^l)$.</figcaption></figure>

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

SWTBLI (error on velocity field): Error on the DA-RANS-p1 velocity field $v_h^{DA}$, wall pressure coefficient $C_{f,h}^{DA}$ and skin friction coefficient $C_{\epsilon,h}^{DA}$ normalized with the respective SA-RANS-p1 physical quantities and the $L^2$-norm of the optimized force vector computed on the solid wall.

<table><thead><tr><th></th><th>||v<sup>DA</sup><sub>h</sub> - v<sub>m</sub>||<sup>2</sup><sub>2:Ω<sub>h</sub></sub><br>in cost function</th><th>||C<sup>DA</sup><sub>f,h</sub> - C<sub>p,x</sub>||<sup>2</sup><sub>2:(Ω<sub>h</sub>∩Ω<sub>p</sub>)</sub><br>a posteriori</th><th>||C<sup>DA</sup><sub>f,h</sub> - C<sub>f,m</sub>||<sup>2</sup><sub>2:(Ω<sub>h</sub>∩Ω<sub>p</sub>)</sub><br>a posteriori</th><th>||f<sub>xyz</sub>||<sup>2</sup><sub>L<sup>2</sup>(Ω<sub>w</sub>)<br>a posteriori</sub></th></tr></thead><tbody><tr><td>L-BFGS</td><td>0.0474</td><td>0.1814</td><td>0.0998</td><td>1.832 × 10<sup>-5</sup></td></tr><tr><td>N-ADAM</td><td>0.0314</td><td>0.7108</td><td>0.0963</td><td>2.072 × 10<sup>-2</sup></td></tr><tr><td>V-N-ADAM-DG</td><td>0.0447</td><td>0.2009</td><td>0.0958</td><td>1.714 × 10<sup>-5</sup></td></tr></tbody></table>

where $\delta_0$ represents the boundary layer thickness at the domain inlet. The RANS equations are coupled with the low Reynolds number Spalart-Allmaras model (referred to as SA-RANS).

### 5.3.1. Mean velocity field reconstruction

In this field-inversion problem, the error metric on the mean velocity fields is adjusted using a geometric weighting function $w(x)$ to emphasize flow features near the wall:

$$ \mathcal{J}(v) = \frac{1}{2} \|w(x)(v_h - v_m)\|_{L^2(\Omega_h)}^2 \quad \text{with} \quad w(x) = \frac{\delta_0}{d_w + 0.1\delta_0} \qquad (37) $$

where $d_w$ is the wall distance used in the SA model. While incorporating the velocity field into the objective function improves velocity field reconstruction, it does not guarantee full mean flow field reconstruction. The density field is obtained by solving the continuity equation, since no corrective terms are applied to this PDE. However, the total energy (or pressure) field remains unconstrained and does not necessarily match the high-fidelity reference, as it is not directly included in the optimization problem. The SA-RANS and DA-RANS mean velocity fields are compared with the DNS reference in Fig. 6, and the optimal forcing components in the momentum equations are presented in Fig. 7 for the three optimizers. In this scenario, L-BFGS and V-N-ADAM-DG show similar results, while the N-ADAM optimum includes a stronger vertical forcing correction. The differing convergence behaviors of the two stochastic optimizers can be attributed to the fundamental differences in their update mechanisms. In the N-ADAM algorithm, the DOFs associated with each component $1 \le r \le N_r$ are updated uniformly, which directs the optimization toward solutions where $f_x$ and $f_y$ share a similar range of values (see Fig. 7c and d). In contrast, the V-N-ADAM-DG algorithm's flexibility enables independent optimization of each $r$th component of the correction term, as updates are guided solely by the gradient of $J_h$ with respect to the $r$th component of the forcing term. Furthermore, the gradient normalization within the V-N-ADAM-DG update rule stabilizes the optimization process, especially in cases where the objective function's sensitivity to one of the $N_r$ components is significantly stronger than to others. Here, a discussion on the regularity of the optimal forcing term is less critical, given that the $p = 1$ discretization provides only $N_p = 3$ modes per component of $f_{DA}$. However, similar behavior to that observed in the L-CYL test case is noted.

Evaluating the quality of the optimal solutions requires the use of an a posteriori error metric on physical quantities not directly included in the objective function, such as the non-dimensional wall pressure and skin friction coefficients. Table 1 provides the mean error metrics, based on the velocity field, pressure $C_p$ and skin friction $C_f$ coefficients, along with the $L^2$-norm of the inferred momentum forcing field $f_{xyz}$ evaluated over the solid wall boundary $\partial\Omega_w$: Although N-ADAM converges to a lower minimum in the objective function space, its error in the wall pressure field is larger. Furthermore, the magnitude of the inferred wall correction is approximately three orders of magnitude larger than that obtained with the other two optimizers, indicating a less physically consistent near-wall behavior.

Visualizations of the wall quantities in Fig. 8a and b clearly reveal a degradation in the pressure field over large sections of the domain, particularly in the boundary layer and recovery zones. It appears that the lower objective function error was achieved by
<figure><img src="image_8.png" /><figcaption>Fig. 6. SWTBLI (error on velocity field): Mean streamwise velocity contour $\frac{u}{u_\omega} \in [-0.1, 0.9]$. Black continuous lines refer to the DNS solution, red dashed lines refer to the baseline solution (a) and blue dashed lines refer to the DA-RANS solution with the L-BFGS (b), N-ADAM (c) and V-N-ADAM-DG (d) optimizers. (a)Baseline; (b)L-BFGS; (c)N-ADAM; (d)V-N-ADAM-DG.</figcaption></figure>

compromising unobserved quantities, resulting in a loss of physical consistency in the optimal forcing term. Ideally, the forcing term $\rho f_{xyz}$ should only account for RFV correction and should not influence solid wall regions, as turbulence is absent there. The vertical profiles of the two components of the corrective forces are shown in Fig. 9 for three stations: the boundary layer region ($x = 25$), the separation region ($x = 30$), and the recovery zone ($x = 35$). At almost every station, L-BFGS and V-N-ADAM-DG exhibit appropriate behavior, with the forcing term approaching zero at the wall. However, the N-ADAM optimum converges to a non-zero value at the wall, indicating that its forcing term not only accounts for turbulence correction but also compensates for errors introduced by pressure gradient mismatches.

### 5.3.2. Mean velocity and wall pressure reconstruction

To better account for the near-wall behavior of the forcing term, an additional term is introduced in the objective function. This new metric specifically targets the mismatch in the wall pressure coefficient:

$$ \mathcal{J}(v) = \frac{1}{2} \|w(x)(v_h - v_m)\|_{L^2(\Omega_h)}^2 + \frac{1}{2} \|C_{p,h} - C_{p,m}\|_{L^2(\partial\Omega_w)}^2 \quad (38) $$

In this field inversion problem, the full mean flow state is constrained at the wall, as the conservative variables can be reconstructed using the velocity and pressure fields, along with the continuity equation. However, the addition of this term increases the complexity of the objective function, making the optimization problem more challenging for the L-BFGS method, which fails to recover the correct mean flow field (the error metrics are larger in Table 2). The V-N-ADAM-DG optimizer appears to be less sensitive to this increase in complexity, likely due to its stochastic nature, which makes it better suited for optimizing complex objective functions. Table 2 also reports the wall-integrated magnitude of the forcing correction, which quantifies deviations from the expected asymptotic decay to zero near the wall. The N-ADAM optimizer exhibits a wall-integrated forcing of approximately $10^{-1}$, indicating pronounced non-physical behavior, whereas both V-N-ADAM-DG and L-BFGS yield values between $10^{-6}$ and $10^{-5}$, consistent with the correct near-wall asymptotic trend.

Fig. 10 compares the near-wall behavior of the forcing correction obtained with the three optimizers. Although N-ADAM benefits from the same stochastic properties and significantly improves the error metric for the wall pressure coefficient, it introduces even
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** SWTBLI (error on velocity field): Streamwise and vertical corrective forces colored contour and stream-wise velocity black isolines for the L-BFGS (a,b), N-ADAM (c,d) and V-N-ADAM-DG (e,f) optimizers. (a) $\rho f_x$ (L-BFGS); (b) $\rho f_y$ (L-BFGS); (c) $\rho f_x$ (N-ADAM); (d) $\rho f_y$ (N-ADAM); (e) $\rho f_x$ (V-N-ADAM-DG); (f) $\rho f_y$ (V-N-ADAM-DG).</figcaption></figure>

stronger forcing terms at the wall. The discrepancy between good minimization results and the expected physical behavior of the forcing term at the wall may be attributed to the non-linear nature of the gradient update. Specifically, the gradient $g$ should ideally approach zero near the wall. However, for the N-ADAM optimizer, it needs to reach machine precision at the wall so that the sign function in Eq. (22) outputs zero. This issue is absent in the V-N-ADAM-DG update rule (see Eq. (25)), where the gradient information is utilized more effectively to maintain the expected wall distribution.

These results demonstrate that lower global error measures, such as those obtained with N-ADAM, can come at the expense of physical realism. In contrast, the V-N-ADAM-DG optimizer achieves a more favorable balance between error reduction and physical consistency. This makes it the preferred choice for applications where maintaining physically interpretable and stable behavior- particularly in near-wall regions of shear flows-is critical for robust predictions of quantities not directly included in the cost function. Another important aspect is that N-ADAM loses the adjoint-consistency properties discussed in Section 3.1 due to the presence of a non-zero forcing term at the wall. In contrast, V-N-ADAM-DG better preserves the adjoint consistency of the problem, ensuring the potential use of the corrected model in adjoint-based CFD applications.
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** SWTBLI (error on velocity field): Skin friction coefficient $C_f$ and dimensionless wall pressure $\frac{\rho_w}{\rho_{w0}}$. Comparison of the reconstructed wall quantities in the case of DA with velocity field data with the SA-RANS (dashed red line) and the DNS (black dots) solutions. DA-RANS quantities results from optimization carried by the L-BFGS (blue long dashed lines), N-ADAM (green dot-dashed lines) and V-N-ADAM-DG (orange double dot-dashed lines) optimizers. (a)Skin friction; (b)Wall pressure.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** SWTBLI (error on velocity field): Vertical profiles of the streamwise and vertical components of the corrective term $\rho f_{xyz}$ at four different x-stations (defined in sub-captions) for the L-BFGS (blue long dashed lines), N-ADAM (green dot-dashed lines) and V-N-ADAM-DG (orange double dot-dashed lines) optimizers. (a)$x_s = 25$; (d)$x_s = 30$; (c)$x_s = 35$; (d)$x_s = 25$; (e)$x_s = 30$; (f)$x_s = 35$.</figcaption></figure>

### 5.4. Shock-wave/turbulent boundary layer interaction (SWTBLI): Turbulent transport equation correction term

We analyze the performance of the three optimizers on the second formulation proposed in Section 4.3 for the objective function (38). The inverse problem is solved for the SWTBLI case, and the optimization results are presented in Table 3. In this formulation, the source term introduced in the SA model aims to correct the eddy viscosity prediction and improve the DA-RANS results while preserving the validity of the Boussinesq hypothesis. The SA model is a turbulence model designed to be robust and well-calibrated
<table><caption>**Table 2**
SWTBLI (error on velocity and wall pressure field): Error on the DA-RANS-p1 velocity field $v_h^{DA}$, wall pressure coefficient $C_{p,h}^{DA}$ and skin friction coefficient $C_{f,h}^{DA}$, normalized with the respective SA-RANS-p1 physical quantities and the $L^2$-norm of the optimized force vector computed on the solid wall.</caption><thead><tr><th></th><th>||v<sup>DA</sup><sub>h</sub> - v<sub>m</sub>||<sup>2</sup><sub>2|Ω<sub>h</sub>|</sub><br>in cost function</th><th>||C<sup>DA</sup><sub>p,h</sub> - C<sub>p,h</sub>||<sup>2</sup><sub>2|Ω<sub>p,h</sub>|</sub><br>in cost function</th><th>||C<sup>DA</sup><sub>f,h</sub> - C<sub>f,h</sub>||<sup>2</sup><sub>2|Ω<sub>f,h</sub>|</sub><br>a posteriori</th><th>||f<sub>xyz</sub>||<sup>2</sup><sub>2|Ω<sub>xyz</sub>|</sub><br>a posteriori</th></tr></thead><tbody><tr><td>L-BFGS</td><td>0.0898</td><td>0.2971</td><td>0.3191</td><td>4.399 × 10<sup>-6</sup></td></tr><tr><td>N-ADAM</td><td>0.0310</td><td>0.0710</td><td>0.0900</td><td>4.227 × 10<sup>-1</sup></td></tr><tr><td>V-N-ADAM-DG</td><td>0.0498</td><td>0.1763</td><td>0.1051</td><td>1.891 × 10<sup>-5</sup></td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** SWTBLI (error on velocity and wall pressure field): Vertical profiles of the streamwise and vertical components of the corrective term $\rho f_{xyz}$ at four different x-stations (defined in sub-captions) for the L-BFGS (blue long dashed lines), N-ADAM (green dot-dashed lines) and V-N-ADAM-DG (orange double dot-dashed lines) optimizers. (a) $x_s = 25$; (b) $x_s = 30$; (c) $x_s = 35$; (d) $x_s = 25$; (e) $x_s = 30$; (f) $x_s = 35$.</figcaption></figure>

<table><caption>**Table 3**
SWTBLI (error on velocity and wall pressure field): Error on the DA-RANS-p1 velocity field $v_h^{DA}$, wall pressure coefficient $C_{p,h}^{DA}$ and skin friction coefficient $C_{f,h}^{DA}$, normalized with the respective SA-RANS-p1 physical quantities. The last column reports the $L^2$-norm of the optimal forcing term computed in the freestream sub-domain $\Omega_{h,\infty} = \{\kappa \in \Omega_h | \chi_S = 1\}$. $\star$: The N-ADAM optimization was terminated due to divergence of the direct RANS solution. The results shown correspond to the last converged DA-RANS iteration prior to divergence.</caption><thead><tr><th></th><th>||v<sup>DA</sup><sub>h</sub> - v<sub>m</sub>||<sup>2</sup><sub>2|Ω<sub>h</sub>|</sub><br>in cost function</th><th>||C<sup>DA</sup><sub>p,h</sub> - C<sub>p,h</sub>||<sup>2</sup><sub>2|Ω<sub>p,h</sub>|</sub><br>in cost function</th><th>||C<sup>DA</sup><sub>f,h</sub> - C<sub>f,h</sub>||<sup>2</sup><sub>2|Ω<sub>f,h</sub>|</sub><br>a posteriori</th><th>f<sub>xyz</sub>(u<sup>DA</sup>) / f<sub>xyz</sub>(u<sup>SA</sup>)<br>total cost function</th><th>||ρf<sub>xyz</sub>||<sup>2</sup><sub>L<sup>2</sup>(Ω<sub>xyz</sub>)</sub><br>a posteriori</th></tr></thead><tbody><tr><td>L-BFGS</td><td>0.2413</td><td>0.1785</td><td>0.8917</td><td>0.2413</td><td>2.610 × 10<sup>-10</sup></td></tr><tr><td>N-ADAM *</td><td>0.1781</td><td>0.2249</td><td>0.3725</td><td>0.1781</td><td>2.170 × 10<sup>-3</sup></td></tr><tr><td>V-N-ADAM-DG</td><td>0.1834</td><td>0.1510</td><td>0.7987</td><td>0.1762</td><td>6.445 × 10<sup>-10</sup></td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>(a) L-BFGS</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) N-ADAM</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) V-N-ADAM-DG</figcaption></figure>

<caption>**Fig. 11.** SWTBLI (error on velocity and wall pressure fields): Source term correction $\rho\tilde{v}f_\nu$ colored contour and stream-wise velocity black isolines for the L-BFGS (a), N-ADAM (b) and V-N-ADAM-DG (c) optimizers. (a)L-BFGS; (b)N-ADAM; (c)V-N-ADAM-DG.</caption>

for boundary layer predictions. However, its performance deteriorates in the presence of wall-separated flows, as observed in this case. The right-hand side of the SA model (A.2) includes a production source term, which is active in shear-stress regions and near the wall, a destruction term to account for turbulence dissipation near the wall, and a diffusion term to ensure smooth predictions at the edge of the turbulence region. In the absence of these three source terms (i.e., in non-turbulent regions such as the far-field), the model ensures that the far-field value of $\mu_t$ remains lower than 3 times the laminar viscosity.

The N-ADAM optimizer is prematurely halted as the direct problem fails to converge after multiple attempts. The L-BFGS optimizer exhibits behavior similar to that observed in the first formulation: the complexity of the objective function, due to the combination of competing error metrics, results in a suboptimal correction. In contrast, the novel stochastic algorithm successfully ensures convergence of the direct problem residuals and achieves a more significant reduction in error metrics, making it the most effective among the three optimization strategies.

As for the previous formulation, the results in Table 3 are complemented by a quantitative metric that aids the interpretation of the inferred correction. Specifically, the $L^2$-norm of the corrective source term $\rho\tilde{v}f_\nu$ is evaluated within the freestream region $\Omega_{h,\infty}$ to quantify deviations from equilibrium in the SA transport equation. The freestream region is defined as $\Omega_{h,\infty} = \{\kappa \in \Omega_h | \chi_S = 1\}$, and is identified through the element-wise sensor proposed in [49]:

$$ \chi_S = \begin{cases} 1, & \text{if } \|S\|_{L^2}^2 < 0.001 \frac{\|v\|_{\infty}\|v_2\|}{\delta_0}, \\ 0, & \text{otherwise,} \end{cases} $$

with $S = \frac{1}{2}(\nabla v + \nabla v^\top) - \frac{1}{d}(\nabla \cdot v)I$ the deviatoric part of the deformation tensor. This metric captures the magnitude of non-physical production or dissipation of eddy viscosity in quiescent flow regions. Among the tested optimization algorithms, N-ADAM yields the largest $L^2$-norm values (on the order of $10^{-3}$), indicating sustained non-equilibrium behavior, whereas both V-N-ADAM-DG and L-BFGS maintain near-equilibrium conditions with significantly smaller magnitudes (on the order of $10^{-10}$).

Fig. 11 compares the distribution of the optimal source term $\rho\tilde{v}f_\nu$ obtained using the three optimizers. The poor performance of the N-ADAM optimizer may again be attributed to its aggressive update rule, which fails to enforce the expected decay of the optimal source term in the far-field region. Hence, the right-hand side terms of the baseline SA model are naturally deactivated, and therefore, the solution in non-turbulent regions should converge to values close to the default far-field eddy viscosity ratio, here set as $\mu_t \le 3\mu$. The inclusion of intrusive corrections in these regions disrupts the equilibrium of the SA model's source terms, leading to spurious production of eddy viscosity in non-turbulent zones. This behavior is confirmed by the eddy viscosity profiles at different
<figure><img src="image_16.png" /><figcaption>**Fig. 12.** SWTBLI (error on velocity and wall pressure fields): Turbulent-to-laminar eddy viscosity ratio profiles at various x-stations (as indicated in the labels) for the L-BFGS (a), N-ADAM (b) and V-N-ADAM-DG (c) optimizers. Solid lines refer to the solution of the DA-RANS simulation and dashed lines refer to the SA-RANS simulation. The legend in the N-ADAM plot (b), identical to those in other figures, has been removed to emphasize the anomalous behavior in the free-stream region. (a)L-BFGS; (b)N-ADAM; (c)V-N-ADAM-DG.</figcaption></figure>

x-stations in Fig. 12, where the DA-RANS results (solid colored lines) are compared with the SA-RANS results (dashed colored lines). The N-ADAM algorithm is the only search algorithm that fails to maintain the expected far-field behavior of the SA model. Moreover, the eddy-laminar viscosity ratio in the far-field reaches nonphysically high values of up to 600, far exceeding the default value of 3 typical of the baseline SA-RANS model (gray dot-dashed isoline).

The optimal distributions obtained using the L-BFGS (Fig. 11a) and V-N-ADAM-DG algorithms (Fig. 11c) are quite similar, except for a higher inferred correction in the recovery zone, downstream of the shock-wave and boundary layer interaction. The slight difference between the two optimal fields suggests that both algorithms are converging toward the same optimal values, but the L-BFGS optimizer terminates earlier due to the complexity of the objective function space. Both optimizers correctly enforce the expected behavior of the turbulence model in non-turbulent regions, as confirmed by the eddy-laminar viscosity ratio profiles at different x-stations in Fig. 12a and c. Both approaches successfully reduce eddy viscosity within the boundary layer region (x = 25) while increasing peak values in the interaction and recovery zones (x ≥ 35).

**6. Conclusion**

In this work, the capabilities and performance of quasi-Newton and stochastic methods have been investigated for physically constrained field-inversion (FI) problems. An analysis of the N-ADAM search algorithm revealed the sharp behavior of the update rule in stochastic optimizers, which can lead to spurious updates of the optimal parameters in regions where the sensitivity of the objective function is relatively low. This makes the use of stochastic algorithms inconsistent with the nature of the inverse problem
they are intended to solve. To mitigate this issue, a novel stochastic optimizer, denoted V-N-ADAM-DG, was developed to ensure smoother updates of the degrees of freedom being optimized, thereby recovering essential properties of the optimization process.

Numerical tests were performed to assess the reconstruction of mean flow fields in compressible RANS simulations using a high-order discretization framework. Two different correction formulations were considered: a vectorial correction applied to the momentum conservation equations, and a scalar source-term correction introduced in the turbulent transport equation. The performance of the L-BFGS, N-ADAM, and V-N-ADAM-DG optimizers was evaluated for several field-inversion (FI) problems in data-augmented (DA)-RANS equations.

The comparison of the resulting optimal forcing terms revealed the limitations of standard stochastic optimizers such as N-ADAM, which exhibited physical inconsistencies in the inferred vectorial corrections near solid walls-regions where a decay to zero is expected. In contrast, the proposed V-N-ADAM-DG optimizer successfully reconstructed the mean flow quantities and ensured that the forcing term converged to zero in regions where the gradient vanishes, such as solid walls and far-field boundaries. A similar behavior was observed for the second correction formulation, where the excessive inferred correction obtained with N-ADAM led to nonphysical production and destruction of the turbulent variable in the free-stream region. This issue was fully mitigated when using the V-N-ADAM-DG optimizer.

A notable advantage of the V-N-ADAM-DG optimizer-absent in quasi-Newton methods like L-BFGS-is its ability to handle highly non-convex objective functions that incorporate multiple physical quantities within the error metric. Furthermore, both the L-BFGS and V-N-ADAM-DG algorithms demonstrated a rapid decay of the DG modal coefficients during polynomial reconstruction due to their specific update rules. This behavior resulted in smoother corrective forcing distributions, a feature not observed in N-ADAM or similar stochastic methods.

The inadequacy of the N-ADAM algorithm for adjoint-based inverse problems with physical constraints was demonstrated both theoretically and through numerical applications. The novel V-N-ADAM-DG algorithm overcomes the shortcomings of N-ADAM and restores the physical consistency of the inverse problem through its improved update formulation. This makes it a promising tool for a broad range of adjoint-based CFD applications, including aerodynamic shape optimization and stability-sensitivity-receptivity analyses.

### CRediT authorship contribution statement

**Bartolomeo Fanizza:** Writing – original draft, Visualization, Validation, Software, Methodology, Investigation, Formal analysis, Conceptualization; **Pedro Stefanin Volpiani:** Conceptualization, Formal analysis, Investigation, Supervision, Writing – review & editing; **Florent Renac:** Conceptualization, Formal analysis, Supervision, Writing – review & editing, Investigation; **Denis Sipp:** Conceptualization, Formal analysis, Supervision, Writing – review & editing, Investigation.

### Data availability

Data will be made available on request.

### Declaration of competing interest

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

### Acknowledgement

This work was partially funded by the French Directorate General for Civil Aviation (DGAC) project LAMA and the internal project MODDA.

## Appendix A. Spalart-Allmaras one-equation model

In this work, we use the Spalart-Allmaras One-Equation model (SA) [31]. The turbulent dynamic viscosity in Section 4.1 is represented as

$$ \mu_t = \rho \tilde{v} f_{v1} \qquad (A.1) $$

with

$$ f_{v1} = \frac{\hat{\chi}^3}{\hat{\chi}^3 + c_{v1}^3}, \quad \hat{\chi} = \frac{\rho \tilde{v}}{\mu}. $$

The turbulent variable $\tilde{v}$ is obtained by the following transport equation

$$ \partial_t \rho \tilde{v} + \nabla \cdot (\rho \tilde{v} \nabla^T) - \frac{1}{\sigma} \nabla \cdot (\xi \nabla \tilde{v} \nabla^T) = c_{b1} \hat{S} \rho \tilde{v} - c_{w1} f_w \frac{\rho \tilde{v}^2}{d_w^2} + \frac{c_{b2} \rho}{\sigma} \nabla \tilde{v} \cdot \nabla \tilde{v}, \qquad (A.2) $$
where the diffusive coefficient reads $\xi = \mu + \rho\tilde{v}$. Moreover, here

$$\hat{S} = S + S^*,$$

with $S = \|\nabla \times v\|_{L^2}$, $S^* = \frac{\tilde{v}f_{v2}}{\kappa^2 d_w^2}$, $d_w$ is the distance to the wall, and

$$f_{v2} = 1 - \frac{\hat{\lambda}}{1 + \hat{\lambda} f_{v1}}, \quad f_w = g \left( \frac{1 + c_{w3}^6}{g^6 + c_{w3}^6} \right)^{\frac{1}{6}}, \quad g = r + c_{w2} (r^6 - r),$$

$$r = \frac{\tilde{v}}{\hat{S}\kappa^2 d_w^2}, \quad c_{w1} = \frac{c_{b1}}{\kappa^2} + \frac{c_{b2}}{\sigma}.$$

For sake of completeness, we give the values of the constants involved in the above expressions: $\sigma = 2/3$, $c_{b1} = 0.1355$, $c_{b2} = 0.622$, $\kappa = 0.41$, $c_{w2} = 0.3$, $c_{w3} = 2$ and $c_{v1} = 7.1$. Negative values of the eddy viscosity (A.1) are limited using a negative version of the SA model [50].

Different modifications of the original SA-model have been presented in order to improve the turbulence closure in a large variety of scenarios. Here, we present the correction proposed for shear flows at low Reynolds numbers [32]. The modified model is identical to the original model, except that the constant $c_{w2}$ is changed to the following function

$$c_{w2,LRe} = c_{w4} + \frac{c_{w5}}{\left[ \left( \frac{\hat{\lambda}}{40} \right) + 1 \right]^2}$$

with $c_{w4} = 0.21$ and $c_{w5} = 1.5$. The effect of $c_{w2,LRe}$ vanishes for large Reynolds numbers value and therefore the calibration of the original model is still preserved for high Reynolds numbers shear flows.

## References

[1] K. Saitou, K. Izui, S. Nishiwaki, P. Papalambros, A survey of structural optimization in mechanical product development, J. Comput. Inf. Sci. Eng. 5 (3) (2005) 214-226. https://doi.org/10.1115/1.2013290

[2] X. Liu, P. Lu, B. Pan, Survey of convex optimization for aerospace applications, Astrodynamics 1 (2017) 23-40. https://doi.org/10.1007/s42064-017-0003-8

[3] J. Martins, A. Ning, Engineering design optimization, Cambridge University Press, 2021. https://doi.org/10.1017/9781108980647

[4] A. Jameson, Aerodynamic design via control theory, J. Sci. Comput. 3 (3) (1988) 233-260. https://doi.org/10.1007/BF01061285

[5] M.B. Giles, N.A. Pierce, An introduction to the adjoint approach to design, Flow, Turbul. Combust. 65 (3) (2000) 393-415. https://doi.org/10.1023/A:1011430410075

[6] J.R. Martins, J.T. Hwang, Review and unification of methods for computing derivatives of multidisciplinary computational models, AIAA J. 51 (11) (2013) 2582-2599. https://doi.org/10.2514/1.J052184

[7] Z. Li, X. Zheng, Review of design optimization methods for turbomachinery aerodynamics, Prog. Aerosp. Sci. 93 (2017) 1-23. https://doi.org/10.1016/j.paerosci.2017.05.003

[8] K. Duraisamy, Z.J. Zhang, A.P. Singh, New approaches in turbulence and transition modeling using data-driven techniques, in: 53rd AIAA Aerospace Sciences Meeting, 2015. https://doi.org/10.2514/6.2015-1284

[9] A.P. Singh, K. Duraisamy, Using field inversion to quantify functional errors in turbulence closures, Phys. Fluids 28 (4) (2016). https://doi.org/10.1063/1.4947045

[10] L. Franceschini, D. Sipp, O. Marquet, Mean-flow data assimilation based on minimal correction of turbulence models: application to turbulent high Reynolds number backward-facing step, Phys. Rev. Fluids 5 (9) (2020). 094603. https://doi.org/10.1103/PhysRevFluids.5.094603

[11] A.S. Cato, P.S. Volpiani, V. Mons, O. Marquet, D. Sipp, Comparison of different data-assimilation approaches to augment RANS turbulence models, Comput. Fluids, 2023. https://doi.org/10.1016/j.compfluid.2023.106054

[12] A. Ferrero, A. Iollo, F. Larocca, Field inversion for data-augmented RANS modelling in turbomachinery flows, Comput. Fluids, 201 (2020). 104474. https://doi.org/10.1016/j.compfluid.2020.104474

[13] B. Fanizza, P. Stefanin, F. Volpiani, E. Renac, D. Martin, Sipp, Adjoint-based optimization for non-linear inverse problems with high-order discretization of the compressible RANS equations, Appl. Math. Model. 142 (2025). 115984. https://doi.org/10.1016/j.apm.2025.115984

[14] C. Zhu, R.H. Byrd, P. Lu, J. Nocedal, Algorithm 778: L-BFGS-B: fortran subroutines for large-scale bound-constrained optimization, ACM Trans. Math. Software 23 (4) (1997). 550-560. https://doi.org/10.1145/279232.279236

[15] A. Mokhtari, A. Ribeiro, Global convergence of online limited memory BFGS, 2014. https://doi.org/10.48550/arXiv.1409.2045

[16] D.P. Kingma, J. Ba, Adam: a method for stochastic optimization, 2014. https://doi.org/10.48550/arXiv.1412.6980

[17] S.Z. Ling, N. Sharp, A. Jacobson, VectorAdam for Rotation Equivariant Geometry Optimization, in: S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, A. Oh (Eds.), Advances in Neural Information Processing Systems, 35. Curran Associates, Inc, 4111-4122, 2022.

[18] F. Bassi, L. Botti, A. Colombo, D.D. Pietro, P. Tesini, On the flexibility of agglomeration based physical space discontinuous Galerkin discretizations, J. Comput. Phys. 231 (1) (2012) 45-65. https://doi.org/10.1016/j.jcp.2011.08.018

[19] A. Björck, Numerical Methods for Least Squares Problems, SIAM, 1996.

[20] L. Franceschini, Modeling strategies for aerodynamic flow reconstruction from partial measurements, Technical Report, Université Paris-Saclay, 2019. Ph.D. thesis, http://www.theses.fr/2019SACLX092.

[21] R. Hartmann, Adjoint consistency analysis of discontinuous Galerkin discretizations, SIAM J. Numer. Anal. 45 (6) (2007) 2671-2696. https://doi.org/10.1137/060665117

[22] T.A. Oliver, L. D, Analysis of dual consistency for discontinuous Galerkin discretizations of source terms, SIAM J. Numer. Anal. 47 (5) (2009) 3507-3525. https://doi.org/10.1137/080721467

[23] R. Hartmann, T. Leicht, Generalized adjoint consistent treatment of wall boundary conditions for compressible flows, J. Comput. Phys. 300 (2015) 754-778. https://doi.org/10.1016/j.jcp.2015.07.042

[24] T. Dozat, Incorporating Nesterov Momentum into Adam, in: Proceedings of the 4th International Conference on Learning Representations, 2016, pp. 1-4.

[25] D.C. Liu, J. On the limited memory BFGS method for large scale optimization, Math. Program. 45 (1-3) (1989) 503-528.

[26] S. J. W. Jorge Nocedal, Quasi-Newton Methods, New York, NY, Springer New York, 2006. https://doi.org/10.1007/978-0-387-40065-5.6

[27] S. J. W. Jorge Nocedal, Line Search Methods, 3 New York, NY, Springer New York, 2006. https://doi.org/10.1007/978-0-387-40065-5.3

[28] J. Sherman, W.J. Morrison, Adjustment of an inverse matrix corresponding to a change in one element of a given matrix, Ann. Math. Stat. 21 (1) (1950) 124-127. https://doi.org/10.1214/aoms/1177729893
[29] C. Xu, J. Zhang, A survey of quasi-Newton equations and quasi-Newton methods for optimization, Ann. Oper. Res. 103 (2001) 213-234. https://doi.org/10.1023/A:1012959223138

[30] S. Wang, Extensions to the Navier-Stokes equations, Phys. Fluids 34 (5) (2022). 053106. https://doi.org/10.1063/5.0087550

[31] P. Spalart, S. Allmaras, A one-equation turbulence model for aerodynamic flows, in: 30th Aerospace Sciences Meeting and Exhibit, 1992, p. 439. https://doi.org/10.2514/6.1992-439

[32] P.R. Spalart, A.V. Garbaruk, Correction to the Spalart-Allmaras turbulence model, providing more accurate skin friction, AIAA J. 58 (5) (2020). 1903-1905. https://doi.org/10.2514/1.J059489

[33] P.S. Volpiani, M. Meyer, L. Franceschini, J. Dandois, F. Renac, E. Martin, O. Marquet, D. Sipp, Machine learning-augmented turbulence modeling for RANS simulations of massively separated flows, Phys. Rev. Fluids 6 (2021). 064607. https://doi.org/10.1103/PhysRevFluids.6.064607

[34] F. Renac, M.D. L.L. Plata, E. Martin, J.B. Chapelier, V. Couaillier, Aghora: A High-Order DG Solver for Turbulent Flow Simulations, IDIHOM: Industrialization of High-Order Methods-A Top-Down Approach: Results of a Collaborative Research Project Funded by the European Union, 2010. https://doi.org/10.1007/978-3-319-12886-3_15

[35] Y. Saad, M.H. Schultz, GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comput. 7 (3) (1986) 856-869. https://doi.org/10.1137/0907058

[36] P.-O. Persson, J. Peraire, Newton-GMRES preconditioning for discontinuous Galerkin discretizations of the Navier-Stokes equations, SIAM J. Sci. Comput. 30 (6) (2008) 2709-2733.

[37] I. Loshchilov, F. Hutter, Decoupled weight decay regularization, 2019. https://doi.org/10.48550/arXiv.1711.05101

[38] D.P. Foures, N. Dovetta, D. Sipp, P.J. Schmid, A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction, J. Fluid Mech. 759 (2014) 404-431. https://doi.org/10.1017/jfm.2014.566

[39] F. Hecht, New development in freefem, J. Numer. Math. 20 (3-4) (2012) 251-265. https://freefem.org/.

[40] H. Wang, S. Xiang, On the convergence rates of Legendre approximation, Math. Comput. 81 (278) (2012) 861-877. https://doi.org/10.1090/S0025-5718-2011-02549-4

[41] J. Déléry, J.-P. Dussauge, Some physical aspects of shock wave/boundary layer interactions, Shock Waves 19 (6) (2009) 453-468. https://doi.org/10.1007/s00193-009-0220-z

[42] D.S. Dolling, Fifty years of shock-wave/boundary-layer interaction research: what next?, AIAA J. 39 (8) (2001) 1517-1531. https://doi.org/10.2514/2.1476

[43] A.J. Smits, J.-P. Dussauge, Turbulent Shear Layers in Supersonic Flow, Springer Science & Business Media, 2006.

[44] G. Settles, D. Dolling, Swept shock/boundary-layer interactions-Tutorial and update, in: 28th Aerospace Sciences Meeting, 1990, p. 375. https://doi.org/10.2514/6.1990-375

[45] J.-P. Dussauge, P. Dupont, J.-F. Debiève, Unsteadiness in shock wave boundary layer interactions with separation, Aerosp. Sci. Technol. 10 (2) (2006) 85-91. https://doi.org/10.1016/j.ast.2005.09.006

[46] S. Pirozzoli, M. Bernardini, Direct numerical simulation database for impinging shock wave/turbulent boundary-layer interaction, AIAA J. 49 (6) (2011) 1307-1312. https://doi.org/10.2514/1.J050901

[47] P.S. Volpiani, M. Bernardini, J. Larsson, Effects of a nonadiabatic wall on supersonic shock/boundary-layer interactions, Phys. Rev. Fluids 3 (8) (2018). https://doi.org/10.1103/PhysRevFluids.3.083401

[48] P.S. Volpiani, M. Bernardini, J. Larsson, Effects of a nonadiabatic wall on hypersonic shock/boundary-layer interactions, Phys. Rev. Fluids 5 (1) (2020). https://doi.org/10.1103/PhysRevFluids.5.014602

[49] W.A.R. Nikhil, J. Oberoi, Larsson, Multi-fidelity parametric sensitivity estimation for large eddy simulation with the Spalart-Allmaras model, J. Turbul. 24 (5) (2023) 195-212. https://doi.org/10.1080/14685248.2023.2212982

[50] S. Allmaras, F. Johnson, P. Spalart, Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model, in: Seventh International Conference on Computational Fluid Dynamics (ICCFD7), 2012, pp. 1-11.