Impact Statement↩︎

This work introduces a fully differentiable framework that integrates data-driven closure models directly into the solution of nonlinear PDE systems. By treating the equations as implicit layers and enabling end-to-end optimization, the approach unifies simulation and learning within a single workflow. This reduces computational duplication and ensures that physical constraints remain embedded throughout training. Applied to the Reynolds-Averaged Navier-Stokes equations for compressible flows, the method demonstrates that differentiable solvers and Neural Networks can be combined in a stable and scalable manner. Beyond this specific case, the framework is broadly applicable to PDE-driven problems with data-informed components, supporting the development of next-generation hybrid modeling strategies at the interface between scientific computing and Machine Learning.

1 Introduction↩︎

Partial Differential Equations (PDEs) are at the heart of mathematical modeling in physics and engineering. They describe how physical quantities such as temperature, pressure, and velocity in fluids evolve in space and time, encapsulating conservation laws such as those of mass, momentum, and energy. The complexity of these equations often renders analytical solutions intractable, especially in realistic geometries or in the presence of non-linearities. As a result, numerical methods have become indispensable tools for solving PDEs, giving rise to a wide range of simulation frameworks across scientific disciplines.

From an engineering perspective, the development of Computational Fluid Dynamics (CFD) codes is crucial for predicting the forces needed to size and manufacture components, especially in high-performance applications such as turbines in jet engines. This process relies on advanced computational and experimental techniques to optimize performance.

The design of aerodynamic profiles is an iterative process, starting with data acquisition, followed by simulations and prototype testing to refine the design. Modern methods integrate optimization algorithms, more recently Machine Learning (ML) [1], [2], to balance aerodynamics, thermodynamics, and structural mechanics, improving the overall design.

In recent years, ML has gained increasing attention as a powerful tool in applied mathematics, particularly for problems traditionally addressed through numerical analysis. Recent works have investigated the integration of ML directly into the numerical discretization of Finite Volume (FV) solvers. Data-driven learned discretizations for hyperbolic conservation laws have been proposed to replace classical gradient or flux reconstructions with Neural Networks (NNs) trained on high-resolution data, achieving improved accuracy on coarse meshes and in the presence of shocks [3], [4]. These approaches effectively modify the internal structure of the CFD solver itself and are typically implemented (including the PDE solver) in dedicated ML-based frameworks. More generally, structure-preserving neural architectures such as volume-preserving transformers have been introduced to ensure stability and long-term consistency in learned dynamical systems [5].

ML’s flexibility in approximating complex functions has led to the development of surrogate models [6]–[8] that can either completely replace or be integrated into and used to accelerate classical solvers, especially when repeated evaluations are required, such as in (shape) optimization or uncertainty quantification [9]. Beyond surrogate modeling, ML techniques have also been used to correct or augment existing numerical models, by learning discrepancies between model predictions and high-fidelity data. In the context of PDEs, NNs have been employed both to represent solutions directly, as in Physics-Informed Neural Networks (PINNs) [10], and to infer hidden physical quantities or closure terms, offering new paradigms for solving inverse problems or discovering governing laws from data [2]. These approaches open new perspectives for hybrid modeling strategies that combine domain knowledge with data-driven learning.

CFD mainly focuses on solving the Navier-Stokes (NS) equations to study the behavior and interactions of fluids. Specifically, for turbulent flows, the Reynolds-Averaged Navier-Stokes (RANS) equations introduce turbulence models to close the system, each model representing a compromise between fidelity and computational feasibility. The choice and accuracy of these turbulence models critically influence the reliability of CFD simulations, especially in the design and optimization of high-performance components.

The accurate simulation of turbulent flows remains one of the major challenges in CFD. Traditional RANS models provide a cost-efficient approach by introducing closure relations for the Reynolds stresses, yet their performance strongly depends on empirical assumptions and model calibration. Classical turbulence closures, such as the \(k\)–\(\varepsilon\), \(k\)–\(\omega\) and Spalart-Allmaras (SA) models [11], while robust and computationally efficient, are known to fail in the prediction of separated or highly anisotropic flows.

The closure modeling problem is a generic feature of nonlinear multiscale PDE systems, arising whenever a reduced or averaged description is sought. A closely related viewpoint emerges in data assimilation, where the full state of a nonlinear multiscale system is reconstructed from partial observations by solving an inverse problem, typically through the introduction of correction or control terms. In variational approaches, these corrections are adjusted to enforce consistency between the model and the data [12]. From this perspective, data assimilation can be interpreted as an implicit closure strategy, in which the inferred corrections encode the effects of unresolved scales or model errors; once learned or parameterized as functions of the resolved state, they can be recast as explicit closure models, thereby bridging inverse problems and reduced-order modeling. [13] provided a comprehensive overview of data assimilation (field-inversion) and ML-augmented RANS modeling, highlighting the opportunities and difficulties in integrating data-driven corrections within physics-based solvers. Furthermore, a more recent view made by [14] emphasized that, despite the promising results, ML-driven turbulence modeling must still tackle issues of generalization, interpretability, and physical consistency.

In the literature, various ML-based approaches have been proposed to enhance RANS closures [15], [16]. For instance, [17] employed NNs to model the effects of variable fluid properties on turbulence transport, while [18] trained data-driven closures for the \(k\)–\(\omega\) SST model, learning discrepancies between RANS and high-fidelity data to improve mean-field predictions. [19] introduced a Deep Learning (DL) model to predict the eddy viscosity field, training directly on the eddy viscosity, yet still achieving improved accuracy in bluff-body simulations, opening the way for turbulent viscosity prediction using NNs. Surrogate modeling has also been applied to turbulence closure, for instance [20] benchmarked surrogate models for RANS closures based on the \(k\)–\(\omega\) SST model.

Despite the progress, most of these methods treat the ML model as an external correction or surrogate, rather than embedding it within an end-to-end differentiable simulation pipeline, as done in [2]. In addition, most ML-RANS studies rely on incompressible solvers. This preference arises because, in the incompressible regime, the Reynolds stress tensor admits an exact decomposition following [21], as commonly exploited in recent works [2], [22]. Furthermore, the incompressible NS equations admit a continuous adjoint formulation, which substantially simplifies the derivation and computation of adjoint-based gradients.

A complementary line of research has explored the use of Graph Neural Networks (GNNs) as data-driven correction models defined directly on unstructured computational meshes. Recent works by Quattromini et al.have proposed physics-constrained GNN architectures for mean-flow reconstruction and data assimilation in CFD [23], [24]. In these methods, the GNN is trained to predict the forcing or closure term that brings a baseline RANS or mean-flow solution into agreement with high-fidelity or experimental data. The graph-based formulation naturally accommodates arbitrary mesh connectivity and allows the use of solver-derived sensitivities during training.

Industrial CFD codes are typically developed over time spans of 10 to 20 years within High-Performance Computing (HPC) environments, often consisting of Python interfaces wrapped around legacy Fortran or C kernels to ensure high performance on large clusters. It is therefore desirable to design a pipeline that enables seamless integration of such CFD solvers with modern ML frameworks like PyTorch, as similarly done by [25].

In the context considered in this work, the PDE solver—specifically a CFD solver—must be treated as a standard PyTorch module, supporting both forward and backward differentiation. At the same time, we aim to preserve the mature, highly optimized numerical algorithms already available in industrial CFD codes—for example, Newton-based nonlinear relaxation schemes and preconditioned GMRes methods [26] for large-scale linear systems. This motivates handling the CFD solver as an implicit layer within PyTorch, with full support for automatic differentiation in both directions.

Furthermore, we aim to introduce ML operators, such as NNs, in an end-to-end structure inside the CFD-PyTorch module. Therefore, once the parameters, or the NN, have been trained, we would still like to rely on the dedicated algorithms embedded in the CFD code to solve the resulting hybrid PDE–ML system efficiently. This requires incorporating the contributions of the trainable model component directly into the Jacobian used for implicit time integration or steady-state solution. The present paper introduces a pipeline that fulfills these objectives.

Here, we demonstrate the approach by coupling PyTorch with the compressible CFD solver BROADCAST [27], which typifies industrial CFD codes built around a Python interface wrapping Fortran kernels. Algorithmic Differentiation (AD) is performed using the software Tapenade [28], enabling efficient evaluation of vector–Jacobian products at a reasonable computational cost [29]. This capability is essential for the practical implementation of adjoint-based optimization.

NNs are adopted as ML models due to the universal approximation theorem [30], which guarantees their ability to approximate a broad class of nonlinear functions. In addition, NNs can be efficiently trained using gradient-based optimization methods. Their capability to handle heterogeneous datasets, together with their flexibility and generality, makes them a robust and widely applicable strategy for processing and predicting data arising from numerical simulations. This approach bridges the gap between physics-based solvers and data-driven learning, extending ML-RANS methodologies towards fully integrated optimization pipelines for aerodynamic design.

This pipeline is applied to the improvement of turbulence models for the RANS equations. In this setting, we introduce an end-to-end differentiable framework that couples the BROADCAST compressible RANS solver with PyTorch, enabling the optimization of NN weights directly within the solver. This approach allows turbulence quantities, such as the eddy viscosity, to be optimized through algorithmic differentiation, thereby providing a fully consistent and automated route for data-driven closure modeling.

The pipeline is first tested on the 2D NASA Wall-Mounted Hump separated flow [31], [32] to optimize a control parameter in the production term of the transport equation of the turbulent variable of the SA model [33], [34]. Subsequently, to further explore the capabilities of the framework, it is applied to the VKI LS-59 turbine blade case [35] with transonic flow to perform an optimization directly on the eddy viscosity, and finally on a dataset generated starting from the benchmark case. This comprehensive approach aims to enhance the predictive accuracy of turbulence models, thereby contributing to more efficient and reliable aerodynamic simulations.

The structure of this paper is as follows. First, the general formulation of the tackled PDEs is introduced, together with a discussion of their abstract representation and closure requirements (§ 2). Next, the proposed optimization framework is described, outlining the possible objective functions and distinguishing between cases where the full state variable is available and cases where it is unknown and must be computed (§ 2.2). The specific problem under investigation is then presented as a benchmark to demonstrate the capabilities of the developed framework, followed by a detailed presentation of the CFD equations in use and closure models (§ 3.1). Subsequently, the technical implementation is discussed, with particular emphasis on the PyTorch-based routines and the design of custom differentiable modules (§ 3.4). The application cases are presented by first describing the considered geometries and then the adopted optimization strategies (§ 4). The generation of the training dataset is described, with particular emphasis on the variations in geometry and flow conditions (§ 4.4.1). Each case is followed by a detailed presentation and analysis of the results obtained with the different approaches. Finally, conclusions are drawn and perspectives for future developments are outlined.

2 Enhancing Partial Differential Equations resolution using Machine Learning↩︎

We start by considering a baseline model \(\mathcal{R}(\boldsymbol{\omega}) = 0\), where \(\boldsymbol{\omega}\) denotes the state variables of interest. Discretization on a mesh using Finite Difference (FD), Finite Element (FE), or FV methods leads to a system \(\boldsymbol{R}(\boldsymbol{w}) = 0\), where \(\boldsymbol{w}\) represents the discretized state. For example, with a FV method with \(m\) variables on a structured mesh of size \(n_i \times n_j\), the total number of degrees of freedom of the state \(N=m \times n_in_j\). To enhance the model accuracy, a discrete additive correction term \(\boldsymbol{f}(\boldsymbol{w}) = \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\), parameterized by coefficients \({\boldsymbol{\vartheta}}\), is introduced. These parameters are determined by minimizing an objective functional \(\mathcal{J}(\boldsymbol{w})\), defined, for instance, from sparse experimental observations or high-fidelity reference data.

Then, the hybrid discrete model can be written as \[\label{eq:generic} \boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) = 0 \, ,\tag{1}\] where \(\boldsymbol{w} \in \mathbb{R}^{N}\), \(\boldsymbol{R}(\boldsymbol{w}): \mathbb{R}^{N} \rightarrow \mathbb{R}^{N}\) and \(\boldsymbol{f}_{{\boldsymbol{\vartheta}}}(\boldsymbol{w}): \mathbb{R}^{N} \rightarrow \mathbb{R}^{N}\). In the following, the implicit time discretization scheme of the hybrid model is assumed to be:

\[\label{eq:timeint} \frac{\boldsymbol{w}_{n+1} - \boldsymbol{w}_{n}}{\Delta t} + \boldsymbol{R}(\boldsymbol{w}_{n+1}) + \boldsymbol{f}_{{\boldsymbol{\vartheta}}}(\boldsymbol{w}_{n+1}) = 0 \, .\tag{2}\]

On the other hand, the correction \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\) may represent either a local or a non-local closure. In the local case, it can be written as \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) = \boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}_{\boldsymbol{\vartheta}} \circ \boldsymbol{\phi}(\boldsymbol{w}))\), where \(\boldsymbol{\phi}(\boldsymbol{w})\) denotes a set of local features and \(\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}\) is a parametric mapping that returns a local correction field \(\boldsymbol{\alpha}\). In the non-local case, the correction takes the form \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) = \boldsymbol{f}\big(\boldsymbol{w}, \boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\big)\), where \(\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}\) depends on the global or neighborhood state and can therefore capture non-local interactions between cells. In practice, \(\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}\) may be instantiated, for instance, as a MultiLayer Perceptron (MLP) in the local setting or as a GNN in the non-local one, although the present formulation remains independent of the specific choice of parametrization.

In the following, we address two algorithmic problems that combine a PDE solver and trainable parameters. First (§ 2.1), for given parameters \({\boldsymbol{\vartheta}}\), we seek the state \(\boldsymbol{w}\) that solves 1 , or: \[\boldsymbol{w} = \mathop{\mathrm{arg\,min}}_{\boldsymbol{w}} \; || \boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) ||^2.\] This requires adapting the solution strategies of the PDE solver underlying \(\boldsymbol{R}\) so that the correction term \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\) is properly taken into account.

Second (§ 2.2), for a given cost functional \(\mathcal{J}(\boldsymbol{w},{\boldsymbol{\vartheta}})\), we seek parameters \({\boldsymbol{\vartheta}}\) such that \[\begin{align} {\boldsymbol{\vartheta}} &= \mathop{\mathrm{arg\,min}}_{\boldsymbol{\vartheta}} \; \mathcal{J}(\boldsymbol{w},{\boldsymbol{\vartheta}}) \\ &\text{s.t. } \boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) = 0. \end{align} \label{eq:min}\tag{3}\] The objective is to optimize the weights \({\boldsymbol{\vartheta}}\) (e.g., within PyTorch), and the key question becomes: what functionalities must the PDE code \(\boldsymbol{R}(\boldsymbol{w})\) provide to enable transparent and smooth forward and backward propagations through the PDE solver?

2.1 Solving for the state for fixed parameters↩︎

For context, it is important to stress that in DL, models are typically constructed as compositions of elementary explicit functions, called layers, where each layer produces its output from its input through a closed-form evaluation. In contrast, the solution of PDEs often relies on the numerical solution of implicit nonlinear systems, where the state is defined as the solution of a residual equation rather than by an explicit formula. This fundamental difference motivates the introduction of implicit operators within differentiable learning frameworks, as done in [36]. Having defined the concepts of explicit and implicit layers, given a fixed set of parameters, the hybrid problem 1 can be solved either by a residual minimization technique (an explicit operation) or a Newton-like relaxation method (an implicit system).

2.1.1 Residual minimization with an explicit layer↩︎

The first strategy to solve problem 1 avoids forming the Jacobian explicitly by recasting the state solve as a residual minimization problem. In this formulation, the state is obtained by minimizing the norm of the residual, typically through an iterative method that requires only the evaluation of the residual itself and back-propagation through the PDE solver \((\partial_{\boldsymbol{w}} \boldsymbol{R})^{\mathsf{T}}\) and through the closure term \((\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}})^{\mathsf{T}}\). Then, the functional to be minimized is fixed as: \[\mathcal{J}(\boldsymbol{w})=\|\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\|^2_Q = \left[\boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right]^{\mathsf{T}}\boldsymbol{Q}\left[\boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right],\] and the differential of \(\mathcal{J}(\boldsymbol{w})\) can be computed as: \[\begin{align} \delta \mathcal{J} &= 2\left[\boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right]^{\mathsf{T}}\boldsymbol{Q}\left[\partial_{\boldsymbol{w}} \boldsymbol{R}(\boldsymbol{w})+\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right] \delta \boldsymbol{w}, \notag \\ &= \left(2\left[(\partial_{\boldsymbol{w}} \boldsymbol{R}(\boldsymbol{w})+\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}))^{\mathsf{T}}\right]\boldsymbol{Q}\left[\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right]\right)^{\mathsf{T}}\boldsymbol{Q}^{-1}\boldsymbol{Q} \delta \boldsymbol{w}, \notag \\ &=\left<2\boldsymbol{Q}^{-1}\left[(\partial_{\boldsymbol{w}} \boldsymbol{R}(\boldsymbol{w}) +\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}))^{\mathsf{T}}\right]\boldsymbol{Q}\left[\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right], \delta \boldsymbol{w} \right>_Q \label{eq:dJ} \end{align}\tag{4}\] where \(\| \cdot \|_Q^2=\left<\cdot,\cdot \right>_Q\) and \(\left<\boldsymbol{w}_1,\boldsymbol{w}_2 \right>_Q=\boldsymbol{w}_1^{\mathsf{T}}\boldsymbol{Q}\boldsymbol{w}_2\) being the discrete inner product, and where \(\boldsymbol{Q}\in\mathbb{R}^{N\times N}\) is a symmetric positive-definite matrix defining this discrete inner-product in the state space. In practice, \(\boldsymbol{Q}\) involves the continuous \(L^2\) inner-product and accounts for the mesh discretization and for the relative scaling of the different state variables. For instance, in a FV discretization, \(\boldsymbol{Q}\) can be chosen as a diagonal matrix containing the cell volumes repeated for each state variable. This definition ensures that the residual norm corresponds to a physically meaningful weighted \(L^2\) norm over the computational domain.

The forward and backward operators can be sketched as: \[\begin{align} {3} \boldsymbol{w} \; &\longrightarrow& \settowidth{\wantedwidth}{\fcolorbox{black}{white}{2\boldsymbol{Q}^{-1}\leftarrow \left(\partial_{\boldsymbol{w}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}} +\partial_{\boldsymbol{w}} \boldsymbol{R}\right)^{\mathsf{T}}\leftarrow \boldsymbol{Q}\leftarrow\left(\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right)}} \makebox[\wantedwidth]{\fcolorbox{black}{white}{ \|\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\|^2_Q }} &\longrightarrow &\mathcal{J} \\ \delta \boldsymbol{w} \; &\longleftarrow& \; \fcolorbox{black}{white}{2\boldsymbol{Q}^{-1}\leftarrow \left(\partial_{\boldsymbol{w}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}} +\partial_{\boldsymbol{w}} \boldsymbol{R}\right)^{\mathsf{T}}\leftarrow \boldsymbol{Q}\leftarrow\left(\boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\right)} \; &\longleftarrow \delta &\mathcal{J}, \end{align}\] where all boxes, handled by a combination of the PDE solver and the automatic differentiation framework (e.g., PyTorch), are explicit since only explicit evaluation operations are involved.

Such an approach, which is similar to PINNs [10], reduces the memory footprint (there is no matrix to invert) and simplifies the implementation, making it attractive when the Jacobian is difficult or expensive to compute. However, it may converge more slowly than Newton-type schemes, and its efficiency depends on the choice of minimization algorithm and the conditioning of the residual landscape.

2.1.2 Relaxation method with an implicit layer↩︎

As an alternative, problem 1 may be solved iteratively by the following relaxation method: \[\label{eq:Newtonstep} \boldsymbol{w}_{n+1}=\boldsymbol{w}_{n}+\delta \boldsymbol{w}, \;\;\;\; \text{with}\;\;\;\; \left(\frac{\mathbf{I}}{\Delta t}+ \partial_{\boldsymbol{w}}\boldsymbol{R}+ \partial_{\boldsymbol{w}}\boldsymbol{f}_{\boldsymbol{\vartheta}} \right)\delta \boldsymbol{w} = -\left(\boldsymbol{R}(\boldsymbol{w}_{n}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_{n})\right).\tag{5}\] The pseudo time-step \(\Delta t\) is small at the beginning of the iterations to accommodate for the initial condition, then it increases to reach \(\Delta t \rightarrow \infty\), where we recover the quadratic convergence properties of the Newton method [37]. The pseudo time-step \(\Delta t\) can be derived from an adaptive stability criterion that evolves with the residual norm, see for example [38] in the case of advection dominated systems with the local \(CFL = c\Delta t/\Delta x\) condition (\(\Delta x\) being the local mesh size and \(c\) the speed of the fastest waves). Note that using a Newton strategy requires the full Jacobian \(\partial_{\boldsymbol{w}} \boldsymbol{R}+ \partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}\) to be invertible at the targeted fixed point, which excludes being in the vicinity of saddle-node, pitchfork or transcritical bifurcation thresholds.

This method requires assembling and inverting the sparse Jacobian by exploiting the stencil-based independence of the numerical scheme. The sparsity pattern of the Jacobian is inferred from the discretization stencil, allowing the use of directional test vectors to probe multiple independent degrees of freedom simultaneously. As a consequence, several Jacobian entries are computed at once using the same test vector, and the corresponding matrix-vector products are performed in parallel, resulting in an efficient sparse assembly procedure.

The Jacobian can be decomposed into two contributions: the part associated with the PDE solver, \(\partial_{\boldsymbol{w}} \boldsymbol{R}\), and the component arising from the linearization of the ML model, \(\partial_{\boldsymbol{w}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\). When \(\boldsymbol{f}_{{\boldsymbol{\vartheta}}}\) is represented by, or contains, a NN, it is essential that the spatial support (or receptive field) of its linear operations does not exceed the stencil width of the underlying numerical scheme. Otherwise, additional couplings would be introduced beyond the assumed sparsity pattern, leading to an incorrect Jacobian structure. This aspect is further discussed in Appendix 9. To relax the constraints on the NN architecture, matrix-free methods could be considered. However, their convergence is typically more challenging without ad hoc preconditioning. This aspect could be explored in future works.

Both methods, explicit and implicit, are compatible with the differentiable framework used in this work, and the choice between them depends on the trade-off between computational cost, implementation complexity, and convergence properties.

2.2 Optimization framework↩︎

In the following, we aim at optimizing the parameters \({\boldsymbol{\vartheta}}\) of a correction term \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\). Similarly to how we previously introduced \(\boldsymbol{Q}\), we introduce a symmetric positive-definite matrix \(\boldsymbol{M}\) defining the natural inner-product in the parameter space: \[\begin{align} \left< {\boldsymbol{\vartheta}}_1, {\boldsymbol{\vartheta}}_2 \right>_M&={\boldsymbol{\vartheta}}_1^{\mathsf{T}}\boldsymbol{M} {\boldsymbol{\vartheta}}_2. \end{align}\]

Two classical situations arise:

1. In the case of data assimilation, we assume \[\label{eq:da} \boldsymbol{f}_{{\boldsymbol{\vartheta}}=\boldsymbol{\alpha}}(\boldsymbol{w})=\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}),\tag{6}\] meaning that the correction \(\boldsymbol{f}_{\boldsymbol{\vartheta}}\) is a function of a volumetric source term \(\boldsymbol{\alpha}\) defined on the computational mesh and handled by the discretization of a continuous operator. In such a case: \({\boldsymbol{\vartheta}}=\boldsymbol{\alpha}\) and the scalar product \(\boldsymbol{M}\) may be chosen to discretize the continuous \(L^2\) inner product, consistently with the inner product used for the state variables.

2. In the case of closure, we assume that \[\label{eq:closure} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})=\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w})),\tag{7}\] meaning that the correction \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\) is a function of the above-introduced parameter \(\boldsymbol{\alpha}\), which now depends on the state \(\boldsymbol{w}\) via a parametric model with parameters \({\boldsymbol{\vartheta}}\), \(\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\). In this case, the natural choice for the scalar product is \(\boldsymbol{M}=\boldsymbol{I}\), since these parameters are not associated with a spatial discretization. Note that: \[\begin{align} \partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})&=\partial_{\boldsymbol{w}}\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w}))+\partial_{\boldsymbol{\alpha}}\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w})) \; \partial_{\boldsymbol{w}} \boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w}), \\ \partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})&=\partial_{\boldsymbol{\alpha}}\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha}_{\boldsymbol{\vartheta}}(\boldsymbol{w})) \; \partial_{\boldsymbol{\vartheta}} \boldsymbol{\alpha}(\boldsymbol{w}), \end{align}\] in which all derivatives are handled by the PDE solver, except for \(\partial_{\boldsymbol{\vartheta}} \boldsymbol{\alpha}(\boldsymbol{w})\), which is treated by the automatic differentiation framework (e.g., PyTorch), as it corresponds to the differentiation of the ML model.

2.2.1 Full-state \(\boldsymbol{w}=\boldsymbol{w}_m\) and explicit layer↩︎

In this section, the full state variable \(\boldsymbol{w}\) is assumed to be known at the training stage, \(\boldsymbol{w}=\boldsymbol{w}_m\). For example, this would be the case if a high-fidelity solver were used to obtain the training data. The optimization problem 3 may then be simplified and replaced by minimizing a cost-functional based on the residue: \[\begin{align} \mathcal{J}({\boldsymbol{\vartheta}}) &:= \| \boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{{\boldsymbol{\vartheta}}}(\boldsymbol{w}_m) \|^2_Q. \end{align}\]

The gradient of the cost-functional can then be obtained straightforwardly: \[\begin{align} \delta\mathcal{J}&=\left<2\boldsymbol{M}^{-1}(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{\boldsymbol{\vartheta}})^{\mathsf{T}}\boldsymbol{Q}\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right), \delta {\boldsymbol{\vartheta}} \right>_{\boldsymbol{M}}. \end{align}\] When using an optimizer based on the \(L^2\) norm, it is essential that the gradient supplied to it lives in a Euclidean space. This inconsistency may result in suboptimal convergence, unstable optimization steps, or biased updates of the parameters. Hence, if \(\boldsymbol{M}\neq I\), we need to introduce the variable \(\tilde{{\boldsymbol{\vartheta}}} = \boldsymbol{N}{\boldsymbol{\vartheta}}\), where \(\boldsymbol{M} = \boldsymbol{N}^{\mathsf{T}}\boldsymbol{N}\) (Cholesky decomposition of a diagonal positive definite matrix), so that \[\begin{align} \left< {\boldsymbol{\vartheta}}_1,{\boldsymbol{\vartheta}}_2\right>_M=\left< \tilde{{\boldsymbol{\vartheta}}_1},\tilde{{\boldsymbol{\vartheta}}_2}\right>_2, \end{align}\] where the subscript \(_2\) refers to the Euclidean norm. The correct gradient to be supplied to the Euclidean optimizer may then be obtained through: \[\begin{align} \delta\mathcal{J} &=\left[2\boldsymbol{M}^{-1}(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{\boldsymbol{\vartheta}})^{\mathsf{T}}\boldsymbol{Q}\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right) \right]^{\mathsf{T}}\boldsymbol{N}^{\mathsf{T}}\boldsymbol{N} \delta {\boldsymbol{\vartheta}} \notag \\ &=\left[2\boldsymbol{N}^{-1}(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{\boldsymbol{\vartheta}})^{\mathsf{T}}\boldsymbol{Q}\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right) \right]^{\mathsf{T}}\delta \tilde{{\boldsymbol{\vartheta}}} \notag \\ &=\left<2\boldsymbol{N}^{-1}(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{\boldsymbol{\vartheta}})^{\mathsf{T}}\boldsymbol{Q}\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right) , \delta \tilde{{\boldsymbol{\vartheta}}}\right>_2. \label{eq:dJN} \end{align}\tag{8}\]

In the present case, if \({\boldsymbol{\vartheta}}\) represents the weights of a neural operator, the forward and backward operators can be sketched as:

\[\begin{align} {3} {\boldsymbol{\vartheta}} &\longrightarrow& \settowidth{\wantedwidth}{\fcolorbox{black}{white}{2\boldsymbol{N}^{-1}\leftarrow \left(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right)^{\mathsf{T}}\leftarrow \boldsymbol{Q}\leftarrow\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right)}} \makebox[\wantedwidth]{\fcolorbox{black}{white}{ \|\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\|^2_Q }} &\longrightarrow &\mathcal{J} \tag{9} \\ \delta {\boldsymbol{\vartheta}} &\longleftarrow& \; \fcolorbox{black}{white}{2\boldsymbol{N}^{-1}\leftarrow \left(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right)^{\mathsf{T}}\leftarrow \boldsymbol{Q}\leftarrow\left(\boldsymbol{R}(\boldsymbol{w}_m)+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_m)\right)} \; &\longleftarrow \delta &\mathcal{J}, \tag{10} \end{align}\] where the boxes are handled by a combination of the PDE solver and the automatic differentiation framework. All layers are explicit in the present case since simple (explicit) evaluation operations are involved in all boxes.

2.2.2 Partial-state and implicit layer↩︎

In the case where the full state variable \(\boldsymbol{w}\) is not known but only partial state information collected in a cost functional \(\mathcal{J}(\boldsymbol{w})\) is available, we can directly tackle the optimization problem in 3 by considering the Lagrangian: \[\label{eq:Lagrangian} \mathcal{L}({\boldsymbol{w}},\boldsymbol{\lambda}, {\boldsymbol{\vartheta}}) = \mathcal{J}({\boldsymbol{w}},{\boldsymbol{\vartheta}}) + \boldsymbol{\lambda}^{\mathsf{T}}(\boldsymbol{R}({\boldsymbol{w}}) + \boldsymbol{f}_{{\boldsymbol{\vartheta}}}(\boldsymbol{w}))\, ,\tag{11}\]

The Karush-Kuhn-Tucker (KKT) optimality conditions are obtained by setting the derivatives with respect to the direct state \(\boldsymbol{w}\) and the adjoint state \(\boldsymbol{\lambda} \in \mathbb{R}^N\) of the Lagrangian to zero: \[\begin{align} \partial_{\boldsymbol{\lambda}} \mathcal{L} \delta \boldsymbol{\lambda} &= 0 ~\Rightarrow~ \boldsymbol{R}(\boldsymbol{w}) + \boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}) = 0 \,, \tag{12} \\ {\partial_{\boldsymbol{w}} \mathcal{L}} \delta{\boldsymbol{w}}&= 0 ~\Rightarrow~ \boldsymbol{\lambda}^{\mathsf{T}}\left( \partial_{\boldsymbol{w}} \boldsymbol{R} + \partial_{\boldsymbol{w}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}} \right) = -\partial_{\boldsymbol{w}} \mathcal{J} \,, \tag{13}\\ \partial_{\boldsymbol{\vartheta}} \mathcal{L} \delta {\boldsymbol{\vartheta}} & = \delta \mathcal{J} ~\Rightarrow~ \left(\partial_{\boldsymbol{\vartheta}} \mathcal{J} + \boldsymbol{\lambda}^{\mathsf{T}}\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right) \delta {\boldsymbol{\vartheta}} = \delta \mathcal{J} \,, \tag{14} \end{align}\] then, having computed the solution \(\boldsymbol{w}\) of 12 , we can use it to compute the Lagrange multiplier from 13 . Finally, by substituting \(\boldsymbol{\lambda}\) in 14 , it is possible to solve for the optimal value of \({\boldsymbol{\vartheta}}\) with a gradient descent strategy, where the gradients are computed as:

\[\label{eq:dL} \delta\mathcal{J}=\left< \boldsymbol{M}^{-1}\left(\left(\partial_{\boldsymbol{\vartheta}} \mathcal{J}\right)^{\mathsf{T}}- \left(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right)^{\mathsf{T}}\left(\partial_{\boldsymbol{w}} \boldsymbol{R}+\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}\right)^{-\mathsf{T}}\left(\partial_{\boldsymbol{w}} \mathcal{J}\right)^{\mathsf{T}}\right),\delta {\boldsymbol{\vartheta}} \right>_M\,,\tag{15}\] where \(\boldsymbol{M}\) is again the natural inner-product depending on the nature of the optimized parameters \({\boldsymbol{\vartheta}}\).

In the present case, the forward and backward operators of the considered implicit layer can be sketched as: \[\begin{align} {5} {\boldsymbol{\vartheta}} &\longrightarrow \settowidth{\wantedwidth}{\fcolorbox{black}{lightgray}{\boldsymbol{M}^{-1}\leftarrow\left(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right)^{\mathsf{T}}\leftarrow \boldsymbol{\lambda} \leftarrow - \left(\partial_{\boldsymbol{w}} \boldsymbol{R}+\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}\right)^{-\mathsf{T}}}} \makebox[\wantedwidth]{\fcolorbox{black}{lightgray}{ \boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})=0, see \S~\ref{sec:sec:solve}}} &\longrightarrow& \settowidth{\wantedwidth}{\; \delta\boldsymbol{w}} \makebox[\wantedwidth]{\boldsymbol{w}} &\longrightarrow \settowidth{\wantedwidth}{ \fcolorbox{black}{white}{ \left(\partial_{\boldsymbol{w}}\mathcal{J}\right)^{\mathsf{T}}}} \makebox[\wantedwidth]{\fcolorbox{black}{white}{ \mathcal{J}(\boldsymbol{w}) }} &\longrightarrow &\mathcal{J} \tag{16} \\ \delta {\boldsymbol{\vartheta}} &\longleftarrow \fcolorbox{black}{lightgray}{\boldsymbol{M}^{-1}\leftarrow\left(\partial_{\boldsymbol{\vartheta}} \boldsymbol{f}_{{\boldsymbol{\vartheta}}}\right)^{\mathsf{T}}\leftarrow \boldsymbol{\lambda} \leftarrow - \left(\partial_{\boldsymbol{w}} \boldsymbol{R}+\partial_{\boldsymbol{w}} \boldsymbol{f}_{\boldsymbol{\vartheta}}\right)^{-\mathsf{T}}} \; &\longleftarrow& \; \delta\boldsymbol{w} &\longleftarrow \fcolorbox{black}{white}{ \left(\partial_{\boldsymbol{w}}\mathcal{J}\right)^{\mathsf{T}}} &\longleftarrow \delta &\mathcal{J}. \tag{17} \end{align}\] The boxes are handled by the PDE solver and the automatic differentiation framework, the ones in light-gray containing implicit operations (involving complex operations such as a Newton solver or a large-scale matrix inverse), the others being completely explicit. We insist on the fact that the light-grey box in the backward operator is implicit: the operator \((\partial_{\boldsymbol{w}}\boldsymbol{R}+\partial_{\boldsymbol{w}}\boldsymbol{f}_{{\boldsymbol{\vartheta}}})^{-\mathsf{T}}\) should not be interpreted as the application of the explicit inverse of the adjoint matrix (which is generally dense and cannot be stored). Instead, it corresponds to solving the linear system involving the adjoint matrix: \[-(\partial_{\boldsymbol{w}}\boldsymbol{R}+\partial_{\boldsymbol{w}}\boldsymbol{f}_{{\boldsymbol{\vartheta}}})^{\mathsf{T}}\boldsymbol{\lambda} = \delta{\boldsymbol{w}},\] which can be done with either direct methods [39], [40] or iterative strategies [33]. The adjoint system inherits the same sparsity pattern as the Jacobian arising from the linearization of the baseline PDE operator and its eigenvalues and condition number remain unchanged with respect to the original system. Hence, from a mathematical point of view, the complexity in solving the adjoint system should be equivalent to solving the direct one.

In the case where \({\boldsymbol{\vartheta}}\) designates a volumetric source term in the PDE solver, one needs, as in the previous section, to move to a Euclidean space by introducing the variable \(\tilde{{\boldsymbol{\vartheta}}}\) (see § 8).

2.3 The end-to-end differentiable learning pipeline↩︎

Thanks to its flexibility, native support for automatic differentiation, and strong integration with scientific computing libraries, PyTorch [41] was selected as the core framework for implementing the computational routines developed in this work. The overall pipeline is organized around a modular structure composed of three fundamental components responsible for (i) enforcing boundary conditions, (ii) evaluating the discrete residual of the governing equations, and (iii) applying external or model-derived forcing terms. This layered design provides a clean interface between high-level Python code and low-level numerical kernels, while ensuring that each stage of the computation remains physically consistent and computationally efficient.

A central requirement of this architecture is full differentiability. To achieve this, each component is implemented by overriding the forward and backward operations, through PyTorch’s torch.autograd.Function interface¹. This mechanism allows complex, optimized numerical routines to be embedded inside a differentiable computational graph, while remaining fully compatible with PyTorch’s automatic differentiation engine.

The computational flow proceeds by sequentially applying the boundary-condition layer to the input state, evaluating the discrete residual of the governing equations, and finally incorporating the forcing contribution. By composing these custom layers into a standard PyTorch nn.Module, the framework automatically manages gradient propagation and the definition of the backward pass, enabling seamless integration with optimization algorithms and ML architectures.

To solve the nonlinear system associated with the discretized equations, the full computation is embedded within an outer layer implementing a Fixed-Point (FP) or Newton-type iteration strategy. This external layer overrides PyTorch’s default backward behavior by solving an associated adjoint system to compute exact gradients with respect to the trainable parameters. This ensures consistent gradient flow across iterative solves and enables end-to-end learning constrained by general nonlinear governing equations.

While the present work focuses on turbulence modeling within the RANS framework, as detailed in the following section, the abstract structure introduced in 1 is in fact representative of a much broader class of physical problems. In reactive flows, for instance, the correction term \(\boldsymbol{f}(\boldsymbol{w})\) may account for source terms arising from chemical reactions. In porous media, it can represent Darcy-type drag effects, while in radiation or multiphase transport, it may encode closure relations for interfacial fluxes or effective material properties. Comparable formulations also emerge in solid mechanics, particularly in the constitutive modeling of complex materials. More generally, this work addresses the modeling of multiscale systems, where similar challenges arise across a wide range of physical contexts. In kinetic theory, such difficulties appear through moment hierarchies derived from the Boltzmann equation [42]. They are also central to many-body statistical physics via the BBGKY hierarchy [43], and to solid mechanics through homogenization techniques and the derivation of effective constitutive laws [44]. From this perspective, the methodology proposed here can be viewed as a general strategy for incorporating modeling closures into PDE-based systems in a systematic and flexible manner. By leveraging differentiable solvers together with adjoint-based optimization, it enables a consistent and scalable coupling between physics-based models and trainable components. As such, it provides a natural framework for integrating ML into scientific computing workflows governed by PDEs, with potential applications well beyond turbulence modeling.

3 Data assimilation and closure with compressible RANS equations↩︎

In this section, we present the PDE problem under investigation and illustrate how the proposed pipeline can be applied to it. First (§3.1), we present the compressible RANS equations, before we introduce a data assimilation and closure problem with two possible tunable parameters: the eddy viscosity \(\mu_t\) (§3.2) and the production term \(\beta\) of the turbulent variable equation (§3.3).

3.1 Compressible RANS equations with Boussinesq hypothesis↩︎

The steady compressible RANS equations can be written as: \[\label{eq:Rans} \begin{cases} \nabla \cdot (\rho \boldsymbol{U}) = 0 \\ \nabla \cdot (\rho \boldsymbol{U} \boldsymbol{U}+p\boldsymbol{I}) - \nabla \cdot \boldsymbol{\tau} - \nabla \cdot \boldsymbol{\tau}^R=0 \\ \nabla \cdot \left[ (\rho E + p) \boldsymbol{U} \right] - \nabla \cdot \left( \boldsymbol{\tau} \cdot \boldsymbol{U} - \boldsymbol{q} \right) - \nabla \cdot \left( \boldsymbol{\tau}^R \cdot \boldsymbol{U} - \boldsymbol{q_t} \right)=0 \end{cases},\,\tag{18}\] where the density \(\rho\), pressure \(p\), viscous stress tensor \(\boldsymbol{\tau}\) and heat flux \(\boldsymbol{q}\) are Reynolds-averaged quantities, while the velocity vector \(\boldsymbol{U} = [U, \, V]\) and the total energy \(E\) are Favre-averaged quantities [45]. Finally, \(\boldsymbol{\tau}^R\) and \(\boldsymbol{q}_t\) denote the Reynolds stress tensor and the turbulent heat flux, respectively.

To close the system, models relating \(\boldsymbol{\tau}^R\) and \(\boldsymbol{q}_t\) to the state variables are required. Most commonly used approaches rely on the Boussinesq eddy viscosity assumption [46], which states that the Reynolds stress tensor \(\boldsymbol{\tau}^R\) is proportional to the deviatoric part of the mean strain-rate tensor \(\boldsymbol{S}\), with proportionality coefficient given by the eddy viscosity \(\mu_t\). This yields: \[\label{eq:Boussinesq} \boldsymbol{\tau}^R = 2 \mu_t \left(\boldsymbol{S} - \frac{1}{3} (\nabla \cdot \boldsymbol{U}) \boldsymbol{I} \right), \;\;\;\; \boldsymbol{q}_t = \frac{\mu_t}{Pr_t} \nabla T,\tag{19}\] where \(Pr_t\) is the (constant) turbulent Prandtl number and \(T\) the Reynolds-averaged temperature. As discussed in [45], a constant turbulent Prandtl number, \(Pr_t = 0.9\), is used in the present simulations and is considered an acceptable baseline choice for the adiabatic configurations investigated here. Then, by incorporating 19 in 18 , the system reads:

\[\label{eq:RansBussi} \underbrace{\left\{\begin{array}{l} \nabla \cdot (\rho \boldsymbol{U}) \\ \nabla \cdot (\rho \boldsymbol{U} \boldsymbol{U}+p\boldsymbol{I}) - \nabla \cdot \boldsymbol{\tau} \\ \nabla \cdot \left[ (\rho E + p) \boldsymbol{U} \right] - \nabla \cdot \left( \boldsymbol{\tau} \cdot \boldsymbol{U} - \boldsymbol{q} \right) \end{array}\right.}_{(*)} \underbrace{\begin{array}{l} \\ -\nabla \cdot \left[2 \mu_t \left(\boldsymbol{S} - (\nabla \cdot \boldsymbol{U}) \boldsymbol{I}/3 \right) \right] \\ -\nabla \cdot \left[ 2 \mu_t \left(\boldsymbol{S}\cdot \boldsymbol{U} - (\nabla \cdot \boldsymbol{U}) \boldsymbol{U}/3 \right) - \mu_t \nabla T / Pr_t \right] \end{array}}_{(**)} \begin{array}{l} =0 \\ =0 \\ =0 \end{array}\tag{20}\]

In the following, after discretizing these equations, we consider two possible choices for the corrected model, characterized by different state variables \(\boldsymbol{w}\), governing equations \(\boldsymbol{R}(\boldsymbol{w})\) and correction terms \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\).

3.2 Data assimilation and closure with \(\mu_t\)↩︎

The first approach consists of defining the baseline model with the state vector \(\boldsymbol{w}=[\rho, \, \rho U, \, \rho V, \, \rho E ]\) and the discrete version of \((*)\) in 20 for \(\boldsymbol{R}(\boldsymbol{w})\). The correction term \(\boldsymbol{f}(\boldsymbol{w},\boldsymbol{\alpha})\) in 6 then represents the left-over part of the equation, the discretization of term \((**)\), with \(\boldsymbol{\alpha}=\mu_t\). The corrected model finally reads: \[\label{eq:NS} \boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}(\boldsymbol{w},{\mu_t}) = 0\tag{21}\] In the data assimilation setting, \(\mu_t\) is treated as a spatially varying tunable parameter, \(\mu_t(\boldsymbol{x})\). In contrast, the closure problem 7 consists of identifying a functional relationship between the eddy viscosity and the state, \(\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w})\). In Section 4.4, this function will be parametrized using a NN with trainable parameters.

3.3 Data assimilation and closure with \(\beta\)↩︎

An alternative approach consists of incorporating a turbulence model directly into \(\boldsymbol{R}\). Here, we consider the SA model, in which the eddy viscosity \(\mu_t\) depends analytically on an additional variable \(\tilde{\nu}\) (i.e., \(\mu_t(\tilde{\nu})\)) governed by [47]: \[\label{eq:turb} \nabla \cdot (\rho \tilde{\nu} \boldsymbol{U})-\nabla\cdot (\sigma^{-1}(\mu+\rho\tilde{\nu})\nabla \tilde{\nu}) =P(\boldsymbol{w})+D(\boldsymbol{w})+C(\boldsymbol{w}).\tag{22}\] Here, \(P(\boldsymbol{w})\), \(D(\boldsymbol{w})\), and \(C(\boldsymbol{w})\) denote the production, destruction, and cross-diffusion terms, respectively [48]. In this framework, the state vector is augmented as \(\boldsymbol{w}=[\rho, \, \rho U, \, \rho V, \, \rho E, \, \rho \tilde{\nu} ]\), and the baseline operator \(\boldsymbol{R}\) now includes the full \((*)+(**)\) equation 20 , and the turbulence variable equation 22 , together with the analytic relation \(\mu_t(\tilde{\nu})\).

The correction term \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w})\) is introduced as an additional source term acting on the turbulence equation of the form \(\beta P(\boldsymbol{w})\) [33], [34], where \(\beta\) is a scalar parameter and \(P(\boldsymbol{w})\) is the production term appearing on the right-hand side of 22 . The corrected model reads: \[\label{eq:NS95SA} \boldsymbol{R}(\boldsymbol{w})+\boldsymbol{f}(\boldsymbol{w},\beta) = 0 \text{, where } \boldsymbol{f}(\boldsymbol{w},\beta) \text{ is the discrete version of } \left( \begin{array}{c} 0 \\ 0 \\ 0 \\ -\beta P(\boldsymbol{w}) \end{array} \right).\tag{23}\]

In the data assimilation framework 6 , \(\beta\) is treated as a spatially varying parameter, \(\beta(\boldsymbol{x})\). The associated closure problem 7 then consists of identifying a functional relation \(\beta_{\boldsymbol{\vartheta}}(\boldsymbol{w})\).

3.4 The end-to-end pipeline applied to BROADCAST↩︎

The BROADCAST solver [27] was chosen as core CFD solver because of the possibility of automatically differentiating [49] its Fortran (F90) routines thanks to the software Tapenade [28]. Indeed, every routine can be linearized to obtain the Jacobian-vector product (JVP) and the vector-Jacobian product (VJP). The (JVP) and (VJP) are used to construct the sparse representations of the Jacobian and of its transpose, making use of the independent stencil defined by the order of the solver. The solution of the linear systems can then be obtained by a sparse LU-inverse.

The BROADCAST solver is a cell-centred FV solver. Boundary conditions are handled via ghost cells. The full state vector therefore consists of components associated with the physical cells \(\boldsymbol{w}_i\) and components associated with the ghost cells \(\boldsymbol{w}_o\). The boundary conditions are enforced by filling the ghost cells (\(\boldsymbol{w}_o=\boldsymbol{B}(\boldsymbol{w}_i)\)) prior to application of the model \(\boldsymbol{R}_i(\boldsymbol{w}_i,\boldsymbol{w}_o)\) and the correction \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_i,\boldsymbol{w}_o)\). The discrete equations, which are equivalent to Eq. 2 , may be written as: \[\boldsymbol{M}\frac{\boldsymbol{w}_{n+1}-\boldsymbol{w}_{n}}{\Delta t} =\boldsymbol{R}(\boldsymbol{w}_{n+1})+\boldsymbol{P}\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_{n+1}),\] where \[\boldsymbol{w}=\left(\begin{array}{c} \boldsymbol{w}_i \\ \boldsymbol{w}_o \end{array}\right), \;\;\;\; \boldsymbol{M}=\left(\begin{array}{cc} \boldsymbol{V} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{0} \end{array}\right), \;\;\;\; \boldsymbol{R}(\boldsymbol{w})=\left(\begin{array}{c} \boldsymbol{R}_i(\boldsymbol{w}_i,\boldsymbol{w}_o) \\ \boldsymbol{w}_o-\boldsymbol{B}(\boldsymbol{w}_i) \end{array}\right), \;\;\;\; \boldsymbol{P}=\left(\begin{array}{c} \boldsymbol{I} \\ \boldsymbol{0} \end{array}\right),\] In this equation, \(\boldsymbol{V}\) is a diagonal matrix representing the volume of each cell.

The current implementation of the BROADCAST solver employs a FV formulation with the high order Flux-Extrapolated MUSCL (FE-MUSCL) reconstruction described in [50]. Moreover, it relies on structured meshes, either single- or multi-block. This choice is primarily driven by the formulation of the numerical schemes and their implementation in the code, where the connectivity between control volumes is implicitly defined through the \((i,j)\) indexing. Such a structure enables efficient stencil operations, predictable sparsity patterns in the Jacobian matrix, and optimized memory access patterns, which results in good computational performance and facilitates the implementation of high-order discretization schemes.

It should be noted, however, that the methods used in this work, applied here to structured meshes using the BROADCAST solver, could also be implemented for unstructured ones. To do so, one would need to provide the JVP and the sparsity pattern of the Jacobian. The strategy remains general and could be applied in the future to different unstructured solvers. The use of structured meshes therefore represents an implementation choice that favors computational efficiency, regular stencil operations, and accurate boundary-layer resolution, at the cost of reduced geometric flexibility and more complex mesh generation for highly intricate geometries.

Finally, within the BROADCAST solver, the governing equations are adimensionalized using freestream quantities as reference scales, i.e., \(\rho_\infty\), \(U_\infty\), and either the static temperature \(T_s\) or the freestream temperature \(T_\infty\).

The computational flow follows the sketch introduced in Section 2.3. The pipeline is structured around a main module delivering \(\boldsymbol{R}(\boldsymbol{w})\), RANS, composed of three internal layers: BC, NS and FORCE, which respectively handle the imposition of boundary conditions (\(\boldsymbol{w}_o=\boldsymbol{B}(\boldsymbol{w}_i)\)), the evaluation of the baseline model (applying \(\boldsymbol{R}(\boldsymbol{w}_i,\boldsymbol{w}_o)\)), and of the correction term (adding \(\boldsymbol{f}_{\boldsymbol{\vartheta}}(\boldsymbol{w}_i,\boldsymbol{w}_o)\)). Each layer serves as an interface to underlying Fortran routines from the BROADCAST solver [27], ensuring that the operations performed at each stage remain physically faithful and computationally efficient.

The custom layers (BC, NS, and FORCE) can be composed into a PyTorch nn.Module: RANS, which automatically handles the chaining of gradients and the definition of the backward pass, enabling seamless integration with optimization algorithms and NN architectures.

Figure 1 provides a sketch of the RANS module, which combines the BC, NS, and FORCE operators.

To solve the nonlinear system arising from the RANS equations with turbulent closure, the entire computation is wrapped in an outer layer that implements a FP iteration strategy based on Newton’s method.

Figure 2 depicts the details of the two implicit layers appearing in the forward 16 and backward 17 operators, highlighting the hierarchical structure of the implicit solver.

Appendix 6 provides a detailed overview of the PyTorch implementation developed in this work, along with the strategy used to interface the Fortran routines with the Python execution flow. The connection with the Deep Equilibrium Models (DEQ) framework [36] is also discussed. The overall workflow is illustrated through several pseudocode listings.

We emphasize that the development of the implicit layer was the most critical part of this work, as defining a robust and efficient implicit solver involves several aspects that must be carefully considered. Throughout this process, multiple strategies were explored to achieve correct gradient incorporation within the routine. Ultimately, the adopted approach required a complete redefinition of the backpropagation mechanism, allowing for a more accurate and reliable computation of gradients, as detailed in Section 2.2.2 and Appendix 8.

In the case of the turbulent closure problem with the eddy viscosity \(\mu_{t}(\boldsymbol{w})\), a Convolutional Neural Network (CNN) [51], [52] was implemented to represent this relation. Its architecture combines convolutions with varying kernel sizes, with the output excluding ghost cells. This design captures multi-scale features while preserving output consistency and ensuring geometric invariance. The kernel sizes of the convolutions have been chosen so that the sparsity pattern of the correction term \(\boldsymbol{f}(\boldsymbol{w},\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w}))\) appearing in 7 matches the one of \(\boldsymbol{R}(\boldsymbol{w})\). This aspect and the proposed CNN are further discussed in Appendix 9.

4 Results↩︎

The routines are first (§ 4.1) validated on the NASA 2D Wall-Mounted Hump for data assimilation of the \(\beta\) parameter (using both explicit and implicit strategies). We then consider a more ambitious configuration and progressively move toward the development of an algebraic turbulence model based on a CNN:

• First (§ 4.2), we verify that the BROADCAST code, combined with the baseline SA model, yields consistent results for the standard VKI LS-59 turbine blade cascade.

• Then (§ 4.3), we demonstrate that, given the SA solution for the VKI LS-59 case, \(\boldsymbol{w}^{SA}=[\rho,\rho \boldsymbol{U}, \rho E]\), we are able to satisfactorily reconstruct the SA eddy viscosity field \(\mu_t^{SA}\) using the explicit residual minimization strategy by tuning \(\boldsymbol{\alpha} = \mu_t\).

• In § 4.4, we first introduce a custom dataset of flow fields obtained with the SA model by varying both flow and geometric parameters. We then construct a CNN-based algebraic turbulence model \(\mu_{t,\boldsymbol{\vartheta}}(\boldsymbol{w})\) that reproduces the SA eddy viscosity \(\mu_t^{SA}\) over this dataset. To this end, the CNN weights \(\boldsymbol{\vartheta}\) are optimized using the explicit strategy by minimizing the residual of a corrected model on the custom dataset.

Building an algebraic turbulence model capable of reproducing the results of a one-equation SA model is an extremely challenging task, even when using a CNN-based formulation (here with a limited kernel width). It is worth recalling that classical algebraic turbulence models, such as the model by [53], are generally well suited for attached turbulent boundary layers, but perform poorly in separated flows and in wake regions behind profiles, where the turbulent viscosity is significantly mispredicted. In such models, non-locality is only partially accounted for through the distance to the closest wall; however, this does not capture streamwise non-local effects, which are particularly important in wake regions. This limitation notably motivated the introduction of transport-equation-based models, in which streamwise non-locality is represented through the advection of turbulence variables. In the present approach, the CNN-based algebraic eddy viscosity model can be regarded as strongly local, since the distance to the closest wall is not included and the kernel width of the CNN is limited to that of the spatial scheme, i.e., a few grid cells in each spatial direction (see § 9 for more details). As a result, the CNN learns local features based on the resolved state \(\boldsymbol{w}=[\rho,\rho \boldsymbol{U}, \rho E]\) (and its derivatives) to predict the eddy viscosity \(\mu_t\).

In the future, learning a CNN-based \(\beta\) turbulence production correction term from high-fidelity Large Eddy Simulations (LES) data appears to be a promising direction, as it naturally captures streamwise non-locality while identifying local features. Increasing the kernel width provides additional flexibility but may violate the stencil locality of the numerical scheme, introducing extra couplings and modifying the sparsity structure of the associated operators.

4.1 2D NASA Wall-Mounted Hump separated flow↩︎

The first case selected for the application of the framework is NASA’s 2D Wall-Mounted Hump (WH) (no plenum), shown in Figure 3. This case is characterized by a separation of the flow from the body with subsequent reattachment and boundary-layer recovery.

It is a well-established benchmark configuration developed at NASA Langley Research Center for the validation of turbulence models in separated flows. It is widely used within the CFD community because it combines a relatively simple two-dimensional geometry with complex physical phenomena, including adverse pressure gradients, boundary-layer separation, recirculation, and reattachment. The case was designed to provide high-qualityreference data, such as LES, against which RANS and hybrid turbulence models can be systematically assessed.

The geometry, as depicted in Figure 3, consists of a two-dimensional channel-like domain in which a smooth hump is mounted on the lower wall. Upstream of the hump, the lower wall is flat, allowing a turbulent boundary layer to develop under nearly equilibrium conditions. As the flow approaches and passes over the hump crest, it experiences a strong adverse pressure gradient that leads to boundary-layer separation downstream of the crest. This separation produces a recirculation bubble whose size and reattachment location are highly sensitive to turbulence modeling assumptions. Further downstream, the flow reattaches and gradually recovers toward a redeveloping turbulent boundary layer.

The canonical operating conditions correspond to a low-speed subsonic flow with a freestream Mach number of approximately \(M \approx 0.1\) and a Reynolds number on the order of \(Re \approx 9.36 \times 10^{5}\), typically based on a reference length associated with the hump geometry. Under these conditions, compressibility effects remain weak. At the inlet, a turbulent boundary layer consistent with the experimental measurements must be prescribed. The lower wall, including the hump, is treated as a no-slip boundary. The outlet is generally modeled as a pressure outlet with a specified static pressure to ensure the correct overall mass flow.

The primary validation quantities for this case include the surface pressure coefficient \(C_p\), the skin-friction coefficient \(C_f\), and mean velocity profiles at selected streamwise stations. The skin-friction distribution is especially important because it clearly identifies the separation and reattachment points.

A distinctive feature of the distributed two-dimensional grids is the non-flat upper boundary, which exhibits a mild contour rather than being straight. Although the hump is physically located only on the lower wall, the shaped upper wall is intentionally introduced to account for wind-tunnel confinement and blockage effects present in the original experiments. In the experimental setup, the test section had finite height and side plates, which influenced the global pressure distribution within the tunnel. A purely two-dimensional simulation with a flat and distant upper boundary would not reproduce these confinement effects accurately, leading to discrepancies in the pressure gradient imposed on the boundary layer over the hump. Since the separation process is driven by the adverse pressure gradient, even small differences in the global pressure field can significantly alter the predicted separation and reattachment behavior. To avoid performing a fully three-dimensional simulation of the entire wind tunnel, the two-dimensional benchmark geometry incorporates a carefully contoured upper wall that effectively mimics the experimental blockage. In this sense, the upper-wall shape acts as a surrogate for three-dimensional tunnel effects, allowing a two-dimensional computation to reproduce more faithfully the pressure distribution and separated-flow characteristics observed in the experiments.

More information on the case and the physical parameters imposed to run the RANS simulation can be found in [31], [32]. Note that this case is part of the ERCOFTAC Database (Classic Collection). It is listed as Case C.83: Wall-mounted two-dimensional hump with oscillatory zero-mass-flux jet or suction through a slot.

4.1.1 Data assimilation with \(\beta\): partial-state measurements and implicit layer↩︎

In this section, we use the corrected model presented in 23 and the implicit optimization framework described in . The forward implicit solution in the red-framed light-grey box is obtained by iteratively solving 5 with a direct-LU method. The backward implicit solution of the adjoint system uses this same inversion technique. For the definition of the cost functional to optimize, we use as a reference solution the LES simulation run by [54]. Here we consider partial-state measurements based solely on the streamwise and wall-normal velocity to define \[\mathcal{J} = \dfrac{1}{2} \norm{\boldsymbol{U}^{LES}-\dfrac{\left(\rho\boldsymbol{U}\right)^{RANS}}{\rho^{RANS}}}^2_Q \, ,\] and the data assimilation process aims at reconstructing the full state \(\boldsymbol{w}=[\rho, \, \rho \boldsymbol{U}, \, \rho E, \, \rho \tilde{\nu} ]\), i.e., in particular the density, the temperature and the turbulence variable which are considered unknown.

In this benchmark test, the routine was initialized with \(\beta\) set to zero throughout the domain. The BROADCAST solver is launched on a mesh of \(n_i \times n_j = 408 \times 108\) cells in the \(x\) and \(y\) directions respectively (\(n = n_i n_j\)).

The optimization algorithm uses the Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) optimizer [55]. The stopping criterion is reached when the variation in the normalized loss function satisfies a specified tolerance threshold (\(10^{-6}\)). Figure 4 shows the behavior of the normalized objective function.

We denote by \(\beta(\boldsymbol{x})~=~\beta^{opt}(\boldsymbol{x})\) the value of the optimized spatial parameter after the final iteration. It is then possible to study the error fields for the \(U\) and \(V\) components of the velocity. Figure 5 and Figure 6 show the local errors on the discretized mesh for the baseline and the optimized models, respectively. It is important to note that not all of the domain is shown in the figures since the LES reference values were available only in the depicted region.

Moreover, the optimization behavior can be assessed by analyzing the final distribution of \(\beta^{opt}(\boldsymbol{x})\) and its impact on the eddy viscosity. In particular, we can compare the eddy viscosity obtained for the baseline model \(\beta(\boldsymbol{x})~=~0\) with that resulting from \(\beta^{opt}(\boldsymbol{x})\). Figure 7 shows three quantities: the initial eddy viscosity computed with the SA turbulence model, the corresponding eddy viscosity \(\mu_t\) obtained using \(\beta^{opt}\), and the optimized correction field \(\beta^{opt}(\boldsymbol{x})\).

We can see that the turbulent viscosity computed with the optimized value of \(\beta\) has been increased in the recirculation bubble. These results demonstrate the potential of the approach and help identify areas for further refinement and optimization. One can also compare different models in terms of the \(C_f\) skin-friction coefficient and the \(C_p\) pressure coefficient. Figure 8 shows that the obtained results are closer to the LES values than those of the baseline SA turbulence model, thereby confirming that the optimization has yielded more accurate predictions.

It is important to note that, for this study, at every iteration performed by the optimizer a FP problem must be solved using a Newton method, after which another linear system must be inverted to compute the adjoint solution, as detailed in Section 2.2.2. Each iteration using the BROADCAST solver takes approximately one minute, depending on the hardware used, and benefits from multiprocessing to construct the sparse Jacobian matrices.

4.1.2 Data assimilation with \(\beta\): full-state measurements and explicit layer↩︎

In order to fully validate the developed framework, the explicit residual minimization strategy is applied to this same case. As mentioned above, such a strategy requires full-state information \(\boldsymbol{w}_m\). We still consider the corrected model presented in 23 but we now apply the explicit optimization framework described in . The full-state \(\boldsymbol{w}_m\) is taken as the one obtained in the last section for the optimized value of \(\beta\). The optimization algorithm is initialized with \(\beta=0\) in the whole domain. Results are obtained using the Stochastic Gradient Descent (SGD) PyTorch optimizer [57], a tolerance threshold of \(10^{-6}\) being fixed for the variation in the normalized loss function.

This strategy is considerably faster than the implicit one, since instead of performing several Newton iterations, only a single evaluation of a Fortran function is needed. Figure 9 shows the behavior of the normalized objective function.

The final values of \(\beta\) are shown in Figure 10.

While the previous approach required solving many linear systems, this one only requires repeated calls to explicit functions of the CFD Fortran code compiled for Python. Therefore, the optimization using this strategy is much faster, performing thousands of iterations within minutes; however, it requires the full state field.

4.2 Transonic flow through the VKI LS-59 turbine blade cascade: Validation of BROADCAST↩︎

We also consider a transonic flow passing between the blades of the VKI LS-59 turbine cascade [35]. The LS-59 geometry is a well-studied linear turbine blade cascade used both in experiments and numerical simulations to analyze aerodynamic losses, boundary-layer behavior, shock-boundary layer interactions, wake turbulence, and transition phenomena [50], [58], [59]. For this case, either a C-mesh or an H-mesh is commonly used in the literature. These mesh types allow for proper discretization of the boundary layer around the airfoil.

The key aerodynamic phenomena to be addressed include:

• Shock-wave formation on the suction side when the flow becomes locally supersonic, and the associated interaction with the boundary layer, possibly leading to separation or shock-boundary layer interaction.

• Wake development downstream of the trailing edge, including turbulence intensity, wake spreading, and mixing losses. These wake effects influence the overall cascade loss and are sensitive to both Mach number and incoming turbulence.

• Influence of varying the outlet isentropic Mach number (subsonic vs transonic) on the flow features above (shock strength, separation, wake behavior) and on the losses.

A mesh is generated to run a RANS simulation on this case. The BROADCAST solver requires structured meshes, as mentioned earlier; therefore, a single-block structured mesh of the inter-blade space is constructed, as shown in Figure 11. The boundary conditions are set as in the transonic case studied in [58], as depicted in Figure 11.

The following parameters were fixed: (i) mesh sizes and geometry: \(n_i~=~300\), \(n_j~=~120\) (\(n~=~n_i n_j\)), geometry from DOE, (ii) scheme: FE-MUSCL of order 7, (iii) turbulence model: SA, (iv) physical parameters: \(Mach_{in}~=~0.275\), \(Mach_{out}\) from DOE, \(\alpha\) (angle of attack) from DOE, \(Re_{in}~=~2.68~\times~10^5\), \(Re_{out}~=~6.8~\times~10^5\), \(T_{ext}~=~300~K\), \(P_{ext}~=~1.01~\times~10^5~Pa\).

By then solving the RANS equations using the relaxation method (see § 2.1.2), it is possible to compare the results with reference or experimental data. To this end, the isentropic Mach number is computed and compared with the values reported in [58]. The results are shown in Figure 12.

4.3 Transonic flow through the VKI LS-59 turbine blade cascade: Data assimilation with \(\mu_t\) with full-state measurements and explicit layer↩︎

We use the corrected model presented in 21 and the explicit optimization framework described in . The considered correction parameter is the eddy viscosity: \[\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w}) ={\boldsymbol{\vartheta}} \, ,\]

As previously explained, if the full velocity field is available, it is in principle possible to use the residual minimization technique to predict turbulent parameters of the flow. In the present case, given the full state, can we correctly predict the eddy viscosity? To check this, a trainable parameter tensor of the size of the mesh is initialized with a constant value throughout the domain and provided as \(\mu_t\) in the routine. The loss function is composed of two terms: first, the norm of the full residual (laminar baseline model plus turbulent forcing) with the state \(\boldsymbol{w}_m\) corresponding to the solution computed in the previous section (see Fig. 12); second, a regularization term defined as \(\gamma \| \mu_t \|^2_M\), where \(\gamma > 0\), to minimize the values of \(\mu_t\) if not needed (the initial condition starts from a non-zero positive \(\mu_t\) value and should tend to zero in the free-stream for instance).

By running the routine for 2000 optimizer iterations using SGD, both the behavior of the loss function and the final parameter field demonstrate the capabilities of this strategy. Figure 13 shows the evolution of the normalized residual norm (fig. a), the final \(\mu_t\) field (fig. b) and the pointwise absolute error of this solution w.r.t the to the SA one, normalized by the maximal value of the SA solution (fig. c).

We can see that, for this case, the optimization framework accurately recovers the eddy viscosity. Since the state variables used are the result of the SA model, the recovered \(\mu_t\) exhibits the same features expected from that turbulence model. Indeed, the retrieved \(\mu_t\) distribution reproduces the expected behavior, correctly capturing the decay within the viscous sublayers near the walls and the destruction in the outer regions of the boundary layers. In the free-stream, the eddy viscosity is advected downstream without destruction, a characteristic behavior of the SA model.The largest relative errors (\(>7\%\)) are localized in the wake region downstream of the blade, where thin convective structures and strong eddy viscosity gradients increase the sensitivity of the prediction. Nevertheless, the overall distribution and magnitude of the eddy viscosity field remain well reproduced.

4.4 Turbine flow modeling with a CNN for \(\mu_t(\boldsymbol{w})\)↩︎

As outlined previously, our objective is to augment RANS simulations with data-driven corrections. To this end, a representative dataset of simulations is required for training the proposed model, with the VKI LS-59 cascade selected as the baseline configuration. In this section, we first present the dataset (§4.4.1), then use the explicit strategy to learn a predictive model for the eddy viscosity based on full-state fields \(\boldsymbol{w}_m\) from the dataset (§4.4.2).

4.4.1 Generation of the dataset↩︎

When using an ML model, special care must be taken, especially if spatial convolutions are required. To balance the requirements imposed by the solver, a single-block structured mesh discretizing the blade passage is considered. For all parameters, the topology and the number of cells in each direction are kept constant. The suction and pressure sides of the blade are split, positioning them in the upper and lower parts of the computational domain, respectively. This solution addresses the challenges of using ML models but requires particular attention during mesh generation to: (i) correctly discretize the boundary layer, (ii) avoid overly stretched cells at the leading and trailing edges of the blade, and (iii) adapt the domain shape to better align the flow with the generalized computational coordinates. To achieve this, a custom elliptic mesh generator has been implemented. Specifically, by varying the control points used to define the camber line of the blade and the angle of rotation, many different meshes can be generated with the same tool.

To illustrate the effectiveness of the mesh generator, Figure 14 presents the collection of blade profiles generated, in comparison with the original blade geometry and its corresponding control points. These variations highlight the flexibility of the parametrization strategy in exploring a wide design space around the reference configuration. The original profile and control points, extracted from [35], are shown for reference.

The resulting dataset, comprising all generated geometries and corresponding simulation data, has been made openly accessible to the research community and is available at [60]. This resource aims to foster reproducibility and further investigations in the field of turbo-machinery optimization and data-driven turbulence modeling. Moreover, this dataset, together with several others structured in the Plaid format, has been used to assess the performance of different surrogate models. A comprehensive analysis of these applications is provided in [7].

To create the dataset, a Design of Experiments (DOE) was generated using the maximum projection algorithm [61]; Table 1 shows the intervals chosen for the parameters², and the remaining parameters were fixed as detailed in Section 4.2.

In Appendix 7 the dataset is analyzed further, showing its completeness and sensitivity with respect to the parameters used in the DOE; see Table 1.

4.4.2 Closure with \(\mu_t\) using a CNN and explicit layer↩︎

We still use the corrected model presented in 21 and the explicit optimization framework described in .

Optimization with respect to an eddy viscosity model depending on the state \(\boldsymbol{w}\) is now performed: \[\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w}) =\mathcal{M}_{{\boldsymbol{\vartheta}}}(\boldsymbol{w})\, ,\] The model \(\mathcal{M}_{{\boldsymbol{\vartheta}}}\) can be any kind of NN. For this study, a CNN is selected since we are dealing with structured meshes. CNNs exploit spatial correlations through convolutional filters, making them particularly effective at capturing local flow features such as gradients, shear layers, and boundary-layer structures. Moreover, compared to fully connected or recurrent architectures, CNNs require fewer parameters for high-dimensional inputs, improving both training efficiency and generalization.

The dataset generated from BROADCAST simulations in the last section has been split into training (80%) and validation (20%) subsets. This division ensures that the model can be optimized effectively while monitoring its generalization performance on unseen data.

Different optimization algorithms available in PyTorch were tested; the L-BFGS optimizer was ultimately selected to produce the reported results, using a mini-batch size of 100. During training, the cost functional \(\mathcal{J}/\mathcal{J}_0\) decreases by approximately two orders of magnitude, reaching \(\mathcal{J}_f/\mathcal{J}_0 = (4.75 \pm 1.39)\times 10^{-3}\) over 10 independent runs (mean \(\pm\) standard deviation). The corresponding validation value is \(\mathcal{J}_f^{\mathrm{val}}/\mathcal{J}_0 = (6.35 \pm 0.83)\times 10^{-3}\), indicating limited sensitivityto the random initialization of the network weights and to the mini-batch training procedure. Although L-BFGS is classically formulated in a deterministic full-batch setting, mini-batch variants have also been investigated in the optimization literature [62]. Although these final values are of the same order of magnitude as those obtained in Fig. 13 (a), the latter corresponds to the direct optimization of \(\mu_t\) for a single configuration, for which nearly perfect convergence could be achieved. In the present case, the optimization is carried out over a diverse set of configurations, so that the network parameters must provide a global compromise rather than an exact optimum for a single case.

Figure 15 provides a visualization of the model predictions for two randomly selected configurations (cases 676 and 679, whose simulation parameters are presented in Table 2) drawn from the validation set. For each configuration, two fields are shown to offer a thorough assessment of the model’s performance. The first field represents the predicted eddy viscosity distribution, denoted as \(\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w})\), across the entire domain. The second field illustrates the absolute pointwise percentage error in the eddy viscosity with respect to the SA eddy viscosity obtained with the BROADCAST solver, \(|\mu_{t,{\boldsymbol{\vartheta}}}(\boldsymbol{w})(\boldsymbol{x}) - \mu_{t}^{SA}(\boldsymbol{x})|/\max_{\boldsymbol{x}}{\mu_{t}^{SA}(\boldsymbol{x})}\).

The results show that the CNN-based algebraic turbulence model, despite its locality constraint, is capable of accurately reconstructing the eddy viscosity in the wake of profiles. This can be attributed to the fact that the CNN does not rely on explicit quantities such as the distance to the closest wall, in contrast to classical algebraic turbulence models. Furthermore, the errors are small in the attached regions of the boundary layers, while they remain non-negligible on the suction side of the profiles near the trailing edge, where flow separation occurs. This clearly indicates that a local algebraic turbulence model does not possess sufficient degrees of freedom to accommodate the full range of flow regimes, including attached boundary layers, separated boundary layers, and wake regions.

To further assess the robustness of the proposed framework, we consider an ensemble of models differing only by random initialization and stochastic optimization. The ensemble is constructed from the same 10 runs used in the training and validation loss analysis, allowing us to quantify prediction variability and estimate epistemic uncertainty.

For case 676, previously presented, we compute the ensemble mean and standard deviation of the predicted turbulent viscosity field. The ensemble mean, shown in Figure 16 (a), closely matches the previously reported prediction, demonstrating the consistency of the model across different training instances.

In contrast, the standard deviation, reported in Figure 16 (b), reveals localized regions of increased variability. Notably, these regions correspond to areas where larger prediction errors were observed in the previous analysis. This correlation between ensemble spread and prediction error suggests that the standard deviation can serve as a meaningful indicator of model uncertainty and reliability.

Overall, this analysis highlights the benefit of ensemble-based approaches for identifying regions of reduced confidence in the predictions, providing additional insight that is not accessible from a single model realization.

Several directions could be explored to further improve the proposed framework. These include the use of richer architectures (e.g., increased channel capacity or alternative field embeddings such as implicit neural representations), the incorporation of tags in the input fields to explicitly encode the boundary conditions information and modifications of the loss function, for instance through importance sampling based on velocity gradients or by adding a regularization term directly on the turbulent viscosity. Finally, much care must be taken to ensure that the trained model behaves robustly when embedded within a Newton solver, in particular with respect to FP stability. These aspects will be addressed in future work.

5 Discussions and conclusions↩︎

The main objective of this work was to demonstrate the capability of the proposed framework to perform end-to-end optimization within a fully differentiable and physically consistent pipeline.

After introducing the optimization problem, the reference configurations, and the training dataset, the first two test cases show that both optimization strategies are fully operational. In particular, the explicit residual minimization and the implicit Newton-like solver correctly propagate gradients through the entire computational chain in a truly end-to-end fashion. This validates the correctness of the proposed formulation and its seamless integration within a differentiable programming environment.

The results further confirm the effectiveness of the proposed approach in predicting the eddy viscosity directly from the conservative flow variables. These findings highlight the potential of the framework as a reliable tool for data-assisted turbulence modeling and provide a solid basis for future developments aimed at improving accuracy, robustness, and scalability.

While demonstrated here for RANS turbulence modeling, the proposed framework is general and extends naturally to PDE-constrained systems with embedded closure models. It thus establishes a consistent and scalable paradigm for integrating trainable components within physics-based solvers via differentiable formulations and adjoint-based optimization.

The most natural extension of this work is the application of the implicit (Deep Equilibrium) strategy to the developed dataset, still targeting the prediction of the eddy viscosity in the Spalart-Allmaras turbulence model. Further extensions may include the adoption of more complex neural architectures, the treatment of different closure problems, and the investigation of transfer capabilities across distinct flow regimes. These directions will contribute to strengthening the synergy between numerical modeling and data-driven methodologies in computational physics.

References↩︎

---

1. It is also possible to define the Jacobian-vector product (jvp) inside a torch.autograd.Function, although it is not needed for optimization purposes; it may, however, be useful to construct the Jacobian for a Newton method.↩︎

2. Further details can be found at [60].↩︎

3. Since the mass matrix \(\boldsymbol{M}\) is diagonal, the associated Cholesky factor \(\boldsymbol{N}=\sqrt{\boldsymbol{M}}\) is diagonal and thus symmetric.↩︎