Addressing a posteriori performance degradation in neural network subgrid stress models

Andy Wu; Sanjiva Lele

doi:10.1017/jfm.2026.11995

Addressing a posteriori performance degradation in neural network subgrid stress models

Published online by Cambridge University Press: 08 September 2026

Andy Wu

and

Sanjiva Lele

Show author details

Andy Wu*: Affiliation:
Department of Aeronautics and Astronautics, Stanford University, Stanford, CA, USA
Sanjiva Lele: Affiliation:
Department of Aeronautics and Astronautics, Stanford University, Stanford, CA, USA Department of Mechanical Engineering, Stanford University, Stanford, CA, USA
*: Corresponding author: Andy Wu, awu1018@stanford.edu

Article contents

Abstract
Introduction
Numerical details
Filtering and neural network details
A priori analysis
A posteriori analysis
Additional a posteriori evaluations
Combined a posteriori analysis
Conclusion
Funding
Declaration of interests
References

Rights & Permissions

Abstract

Content of image described in text.

Neural network subgrid stress models often have a priori performance that is far better than the a posteriori performance, leading to neural network models that look very promising a priori completely failing in a posteriori large eddy simulations (LESs). This performance gap can be decreased by combining two different methods, training data augmentation and reducing input complexity to the neural network. Augmenting the training data with two different filter shapes before training the neural networks has no performance degradation a priori as compared with a neural network trained with one filter. A posteriori, neural networks trained with two different filters are more robust across two LES codes with different high-order numerical schemes. In addition, by ablating away the higher-order terms input into the neural network, the a priori versus a posteriori performance changes become less apparent. When combined, neural networks that use both training data augmentation and a less complex set of inputs have a posteriori performance more reflective of their a priori evaluation.

JFM classification

Turbulent Flows: Turbulence modelling Mathematical Foundations: Machine learning

Information

Type: JFM Papers
Information: Journal of Fluid Mechanics , Volume 1042 , 10 September 2026 , A25

DOI: https://doi.org/10.1017/jfm.2026.11995 [Opens in a new window]
Copyright: © The Author(s), 2026. Published by Cambridge University Press

1. Introduction

Large eddy simulation (LES) provides an economical paradigm for high-fidelity flow prediction by resolving the larger, energy containing scales and modelling the spatial scales of turbulence that are smaller than the grid scale through a subgrid-scale (SGS) closure (Sagaut Reference Sagaut2006). Spatially filtering the Navier–Stokes equations results in an unclosed term called the subgrid stress tensor, $\tau _{\textit{ij}}$ , which represents the effect of the subgrid spatial scales of turbulence on the resolved flow (Pope Reference Pope2000).

Recently, with the availability of direct numerical simulation (DNS) data online and advances in deep learning theory, neural networks have been explored as subgrid stress models (Sarghini et al. Reference Sarghini, de Felice and Santini2003; Beck et al. Reference Beck, Flad and Munz2019; Xie et al. Reference Xie, Wang and Weinan2020a ,Reference Xie, Yuan and Wang b ; Park & Choi Reference Park and Choi2021; Stoffer et al. Reference Stoffer, van Leeuwen, Podareanu, Codreanu, Veerman, Janssens, Hartogensis and van Heerwaarden2021; Cheng et al. Reference Cheng2022; Kang et al. Reference Kang, Jeon and You2023; Cho et al. Reference Cho, Park and Choi2024; Maejimam & Kawai Reference Maejimam and Kawai2024; Wu & Lele Reference Wu and Lele2025b ). A persistent challenge is the discrepancy in neural network a priori and a posteriori performance (Beck et al. Reference Beck, Flad and Munz2019; Beck & Kurz Reference Beck and Kurz2021; Park & Choi Reference Park and Choi2021; Stoffer et al. Reference Stoffer, van Leeuwen, Podareanu, Codreanu, Veerman, Janssens, Hartogensis and van Heerwaarden2021), which can be interpreted as a distribution shift between the training data and the LES simulation (Hu et al. Reference Hu, Subramaniam, Kuang, Lin, Yu, Hannah, Brenowitz, Romero and Pritchard2025). Due to this distribution shift, ad hoc methods are sometimes used, such as adding an additional Smagorinsky term to help the neural network dissipate energy so that simulations remain stable (Beck et al. Reference Beck, Flad and Munz2019). Two mechanisms may contribute to the distribution shift. First, the input fields provided to the neural network during training, which are generated assuming an explicit filter (oftentimes a box filter), can differ from the LES resolved fields, since the LES filter is implicit and is a product of both the grid and the numerical method (Sagaut Reference Sagaut2006). Second, the output SGS stress statistics also differ across filtering operations. For example, the omission or use of two-thirds dealiasing in spectral methods modifies the effective filter transfer function, producing different SGS statistics. This motivates the hypothesis that neural networks trained on one filter may struggle to generalise to LES solvers a posteriori. To address these two mechanisms of distribution shift, this work introduces a multi-filter data augmentation strategy that exposes the neural network to various plausible filtered inputs and SGS stress distributions. These sources of distribution shift are shifts in the spatial distributions. There are also potential distribution shifts as the neural network subgrid stress model is time advanced in the simulation itself, which can manifest as a source of error accumulation that causes a distribution shift (Sirignano et al. Reference Sirignano, MacArt and Freund2020). This work focuses on addressing the spatial distribution shifts instead of those due to time advancement.

Reinforcement learning approaches (Bae & Koumoutsakos Reference Bae and Koumoutsakos2022; Kim et al. Reference Kim, Kim, Kim and Lee2022; Guan et al. Reference Guan, Subel, Chattopadhyay and Hassanzadeh2023; Kurz et al. Reference Kurz, Offenhäuser and Beck2023) aim to bypass these issues by learning the SGS model (or wall model) directly in a LES solver. However, models learned in this manner are limited in their transferability to other solvers, as they learn the combined effect of solver-specific numerical effects and turbulence physics. Similarly, end-to-end learning approaches, where a differentiable solver is used alongside a neural network model (Duraisamy Reference Duraisamy2021; Sirignano & MacArt Reference Sirignano and MacArt2023; Kim et al. Reference Kim, Shankar, Viswanathan and Maulik2024) can suffer similar issues, as the neural network learned in this setting is also unlikely to transfer directly to different flow solvers, having been optimised for a specific flow solver and obtaining gradients that are impacted by the flow solver’s numerical method. In addition, the training cost of the neural network is much higher, as it requires LES to be run during training. In contrast, the present supervised learning approach aims to not learn solver peculiarities, and incorporates known filtering operations into the training data. The filters used do not mimic the specific transfer function of any specific LES solver and instead are inspired by common LES numerical behaviour. However, many works that do adopt a supervised approach where the filtered DNS data are the ground truth will also not learn solver peculiarities. Instead, these works only use one filter (oftentimes a box filter) when assembling the training data. Therefore, instead of learning solver peculiarities, the neural networks can now overfit to filter peculiarities, causing a distribution shift a priori versus a posteriori.

Furthermore, another potential source of distribution shift can come from the use of increasingly complex neural network inputs, such as higher-order powers of the velocity gradient tensor ( $V_{\textit{ij}}$ ). For example, the second invariant $Q = 0.5(P^2-\operatorname {trace}(V_{\textit{ij}}^2)))$ , where $P = \operatorname {trace}(V_{\textit{ij}})$ and $V_{\textit{ij}}$ is the velocity gradient tensor, and the third invariant $R = det(V_{\textit{ij}})$ are more sensitive to aliasing and differences in numerical methods (where $det$ denotes the determinant). Therefore, a systematic ablation of neural network input complexity is performed to analyse the effect of inputs to neural network performance a posteriori.

The contributions of this work are twofold. First, multi-filter (two and three filter) data augmentation is shown to improve the robustness of learned SGS models across small architecture perturbations and across LES solvers with different high-order numerical methods. Second, an ablation of input feature complexity reveals that more complex neural network inputs suffer more from a posteriori performance degradation. Together, these results enable the training of neural networks that are more robust a posteriori.

2. Numerical details

2.1. Large eddy simulation

The LES governing equations involve the filtered Navier–Stokes equations. The filtering operation is defined as

(2.1)

\begin{equation} \bar {\chi }(\boldsymbol{x}, t) = \int \! {G(\boldsymbol{r})\chi (\boldsymbol{x - r}, t) {\rm d}\boldsymbol{r}}, \int \! {G(\boldsymbol{r}) {\rm d}\boldsymbol{r}} = 1, \end{equation}

where $\chi$ is the flow variable to be filtered, $G(\boldsymbol{r})$ is the filtering kernel and $\bar {\chi }$ is the filtered flow variable. In all sections, the overbar operator is used to denote filtered quantities. For robust subgrid stress modelling, the filtering operation must not be oscillatory in physical space, as this would create additional small scales purely through the act of filtering (Sagaut Reference Sagaut2006). This means that the spectral cutoff filter, which is in the shape of a sinc function in physical space (while being sharp in spectral space) is not appropriate as a filter for subgrid stress modelling for general application. With the definition of the filter above, the LES governing equations are written as (in Einstein’s summation notation)

(2.2)

StartLayout 1st Row 1st Column StartFraction partial differential u overbar Subscript i Baseline Over partial differential x Subscript i Baseline EndFraction 2nd Column equals 0 comma 2nd Row 1st Column StartFraction partial differential u overbar Subscript i Baseline Over partial differential t EndFraction plus u overbar Subscript j Baseline StartFraction partial differential u overbar Subscript i Baseline Over partial differential x Subscript j Baseline EndFraction 2nd Column equals minus StartFraction 1 Over rho EndFraction StartFraction partial differential p overbar Over partial differential x Subscript i Baseline EndFraction plus nu StartFraction partial differential squared u overbar Subscript i Baseline Over partial differential x Subscript j Baseline partial differential x Subscript j Baseline EndFraction minus StartFraction partial differential tau Subscript italic ij Baseline Over partial differential x Subscript j Baseline EndFraction comma 3rd Row 1st Column tau Subscript italic ij 2nd Column equals ModifyingAbove u Subscript i Baseline u Subscript j Baseline With bar minus u overbar Subscript i Baseline u overbar Subscript j Baseline comma EndLayout

\begin{equation} \begin{aligned} \frac {\partial \bar {u}_i}{\partial x_i} &= 0 ,\\ \frac {\partial \bar {u}_i}{\partial t} + \bar {u}_{\!j} \frac {\partial \bar {u}_i}{\partial x_{\!j}} &= - \frac {1}{\rho } \frac {\partial \bar {p}}{\partial x_i} + \nu \frac {\partial ^2 \bar {u}_i}{\partial x_{\!j} \partial x_{\!j}} - \frac {\partial \tau _{\textit{ij}}}{\partial x_{\!j}} ,\\ \tau _{\textit{ij}} &= \overline {u_i u_{\!j}} - \bar {u}_i\bar {u}_{\!j} ,\end{aligned} \end{equation}

where $\overline {p}$ is the filtered pressure, $\rho$ is the density, $\overline {u}$ is the filtered velocity and $\nu$ is the kinematic viscosity. As seen, an unclosed term, $\tau _{\textit{ij}}$ , is obtained by filtering the momentum equation. In LES, only the filtered values of the dependent variables are available, thus, $\tau _{\textit{ij}}$ cannot be computed and has to be modelled. In the present work, the neural network will learn to approximate $\tau _{\textit{ij}}$ through filtered DNS data.

2.2. A posteriori flow solvers

Since this work involves drawing conclusions about a posteriori performance, two in-house LES flow solvers, PadeOps and PadeLibs (with different numerical methods) are used (Ghate & Lele Reference Ghate and Lele2017; Song et al. Reference Song, Ghate, Matsuno, West, Subramaniam and Lele2024). PadeOps is a Fourier-spectral code with Runge Kutta (RK)4–5 or RK3 time stepping (Bogacki & Shampine Reference Bogacki and Shampine1996). Therefore, two-thirds dealiasing is applied due to the nonlinear component in the convective acceleration term. Meanwhile, PadeLibs is a sixth-order compact code with RK3 time stepping (with no 2/3 dealiasing). These two codes have different numerical methods, so the a posteriori conclusions obtained are generalisable across these high-order discretisations. As seen in a later section, these two flow solvers create two different environments for the neural network when considering subgrid stress modelling.

3. Filtering and neural network details

Neural network SGS models trained outside of the LES solver and tested a posteriori depend on explicitly filtered DNS data, making the choice of filters critical. Mismatches between the explicit filters chosen and the implicit filters encountered in a posteriori LES can produce both input-side and output-side distribution shifts. To reduce such mismatches, the present work adopts multi-filter training, selecting explicit filters to diversify the training distribution so that the filters span representative LES filtering behaviours without approximating a specific LES solver’s transfer function. The first explicit filter chosen is the box filter, as it is a common filter used in previous studies. The box filter is defined as follows:

(3.1)

\begin{equation} \bar {\chi }(\boldsymbol{x}, t) = \int _V \frac {\chi (\boldsymbol{x - r}, t)}{|V|} \, {\rm d}\boldsymbol{r}, \end{equation}

where $V$ is the volume (or region) of the box filter and $|V|$ is its measure. The box filter is not oscillatory in physical space, but is oscillatory in spectral space. The one-dimensional kernel of the box filter, $G_{\textit{box}}(r)$ , is also given as

(3.2)

upper G Subscript italic box Baseline left parenthesis r right parenthesis equals StartLayout Enlarged left brace 1st Row 1st Column StartFraction 1 Over upper Delta EndFraction 2nd Column for StartAbsoluteValue r EndAbsoluteValue less than or equals StartFraction upper Delta Over 2 EndFraction comma 2nd Row 1st Column 0 2nd Column for StartAbsoluteValue r EndAbsoluteValue greater than StartFraction upper Delta Over 2 EndFraction comma EndLayout

\begin{equation} G_{\textit{box}}(r) = \begin{cases} \frac {1}{\varDelta } & \text{for } |r| \le \frac {\varDelta }{2} ,\\[9pt] 0 & \text{for } |r| \gt \frac {\varDelta }{2} ,\end{cases} \end{equation}

where $\varDelta$ is the filter width.

3.1. Discrete spatial cosine filter: a localised low-pass filter

To obtain a ‘cleaner’ grid cutoff at the Nyquist wavenumber, a discrete spatial cosine filter (DSCF) is used. The DSCF attenuates frequencies faster than the box filter after the grid cutoff while avoiding the oscillatory kernel of an ideal spectral cutoff filter in physical space. The one-dimensional kernel for the localised low-pass filter is given in (3.3) (the authors are grateful to an anonymous referee for suggesting the physical space formulation of the DSCF adopted in the paper. The original treatment in spectral space was equivalent but cumbersome)

(3.3)

upper G Subscript italic DSCF Baseline left parenthesis r right parenthesis equals StartLayout Enlarged left brace 1st Row 1st Column StartFraction 1 Over upper Delta EndFraction left parenthesis 1 plus cosine left parenthesis StartFraction 2 pi r Over upper Delta EndFraction right parenthesis right parenthesis 2nd Column for StartAbsoluteValue r EndAbsoluteValue less than or equals StartFraction upper Delta Over 2 EndFraction comma 2nd Row 1st Column 0 2nd Column for StartAbsoluteValue r EndAbsoluteValue greater than StartFraction upper Delta Over 2 EndFraction period EndLayout

\begin{equation} G_{\textit{DSCF}}(r) = \begin{cases} \frac {1}{\varDelta } \left ( 1 + \cos \left ( \frac {2\pi r}{\varDelta } \right ) \right ) & \text{for } |r| \le \frac {\varDelta }{2} ,\\[9pt] 0 & \text{for } |r| \gt \frac {\varDelta }{2} .\end{cases} \end{equation}

As seen, this filter is a single mode approximation to the spectral cutoff filter (plus the zero-wavenumber component). For a specific filter width, the DSCF cutoff frequency $f_c$ can be found to be $f_c = {1}/{\varDelta }$ , with $\varDelta$ being the filter width. In general, when preparing training data, the cutoff frequency is set to be the grid-scale cutoff for LES. For example, if one wants to filter the DNS data at $\varDelta _{\textit{LES}} = 16 \varDelta _{\textit{DNS}}$ , then the desired $f_{\textit{cd}} = {0.5}/{16} = 0.03125$ (0.5 due to Nyquist frequency), while when filtering at $\varDelta _{\textit{LES}} = 32 \varDelta _{\textit{DNS}}$ , the desired $f_{\textit{cd}} = {0.5}/{32} = 0.0156$ is set. This means that when filtering at $\varDelta _{\textit{LES}} = 16 \varDelta _{\textit{DNS}}$ , the desired $f_{\textit{cd}}$ implies that a 33 point DSCF is needed (which has a filter width of 16 points to the left and 16 points to the right to obtain $f_{\textit{cd}} = f_c$ ), and when filtering at $\varDelta _{\textit{LES}} = 32 \varDelta _{\textit{DNS}}$ , a 65 point DSCF is needed. Both the physical space and spectral space characterisations of the DSCF are shown in figure 1.

Figure 1.

Characterisations of the DSCF in both physical and spectral space. $\omega _c = 2 \pi f_c$ corresponds to the filter cutoff, here designated at a $16\varDelta _{\textit{DNS}}$ filter width, while $\omega = 2 \pi f_n, f_n \in [-0.5, 0.5]$ denotes the normalised angular frequency. Two different DSCF widths are shown, one with a 17 and the other with a 33 point support. The DSCF with the smaller filter width has twice the cutoff frequency.

Two line graphs compare low-pass filter kernels and frequency responses of different filters.

Figure 1. Long description

Panel A: A line graph titled ‘Low-pass filter kernel’ shows the amplitude on the vertical axis and x on the horizontal axis. Two lines represent different filter kernels: DSCF (33-pt) in blue and DSCF (17-pt) in orange. The blue line has a broader peak around zero, while the orange line has a narrower peak. Panel B: A line graph titled ‘Frequency responses of low-pass filters’ displays normalized magnitude on the vertical axis and normalized angular frequency on the horizontal axis. Four lines represent different filters: DSCF (33-pt) in blue, DSCF (17-pt) in dashed orange, Gaussian (16x) in green with dots, and Box (16x) in grey with dots. The blue and dashed orange lines peak at different frequencies, with the blue line having a lower cutoff frequency. The green and grey lines show different frequency responses compared to the blue and orange lines.

As seen in figure 1, the frequency response is not perfectly rectangular, instead exhibiting a smoother roll-off with mild high-wavenumber oscillations when compared with a box filter. Meanwhile, in physical space, the filter is strictly positive, non-oscillatory and has finite support. As seen, the 17 point DSCF has the same filter shape as the 33 point DSCF, but has around $2 \times$ less points in the spatial kernel, and thus has twice the cutoff frequency. The 33 point DSCF matches the initial decay of the box filter but significantly suppresses the oscillations (peak oscillation magnitude of 0.028 vs 0.22). In addition, it attenuates frequencies immediately beyond the grid cutoff more than the Gaussian filter. This is especially seen in the normalised angular frequency range of 0.2–0.4, where the 33 point DSCF is over 4–5 times lower than the Gaussian.

However, one advantage of the Gaussian filter kernel is that it decays higher wavenumbers exponentially, while the DSCF has mild oscillations. While a wider Gaussian filter kernel (and potentially tuning the standard deviation) can match the passband attenuation of the DSCF, near the grid cutoff, the DSCF will still attenuate frequencies faster. This is especially important as filters attenuating frequencies slightly larger than the grid cutoff more aggressively result in less grid-scale noise contaminating the filtered fields. This ensures that the neural network is fed a cleaner signal, which is more free from the grid-scale noise contamination that can hinder convergence in deep learning models.

3.1.1 Two-thirds dealiasing as an additional filter

PadeOps, the Fourier-spectral code with two-thirds dealiasing, uses a filtering operation that modifies the effective transfer function due to the two-thirds dealiasing operation. To study how multi-filter training interacts with this procedure, a filter is derived that is consistent with two-thirds dealiasing. The total subgrid stress can be calculated by applying the two-thirds dealiasing operation on top of the already filtered Navier–Stokes equations, defined as the tilde operation

(3.4)

\begin{equation} \begin{aligned} \tau _{\textit{ij}}^f &= \widetilde {\overline {u_i u_{\!j}}} - \tilde {\bar {u}}_i\tilde {\bar {u}}_{\!j}. \end{aligned} \end{equation}

As seen in (3.4), the subgrid stress tensor, after accounting for this two-thirds dealiasing operation, is denoted as $\tau _{\textit{ij}}^f$ . The two-thirds dealiasing operation is computed as an additional spectral cutoff filter on top of the box filtered data. This filtering operation will be called box and two-thirds (BTF), and is used alongside the box filter to form the augmented dataset for the Fourier-spectral code. Note that the BTF is not intended to exactly replicate the transfer function of Fourier-spectral codes, and it is unlikely that the transfer function of a Fourier-spectral code with 2/3 dealiasing is identical to BTF. Rather, BTF provides the neural network with exposure to the type of behaviour that two-thirds dealiasing qualitatively introduces, to reduce the training data distribution shift a posteriori.

3.1.2. Neural network details

The neural network used is a hybrid tensor basis neural network (TBNN) from Ling et al. (Reference Ling, Kurzawski and Templeton2016) and graph neural network (GNN). In a classic TBNN, an artificial neural network (ANN) maps inputs to scalar coefficients of the tensor basis expansion. Here, the ANN front-end architecture is substituted with a GNN, which allows for learned input stencils instead of prescribed input stencils (Abekawa et al. Reference Abekawa, Minamoto, Osawa, Shimamoto and Tanahashi2023). To recap, the tensor basis expansion for LES, where the neural network predicts 8 scalar coefficients $c_1-c_8$ , can be written as (Stallcup & Dahm Reference Stallcup and Dahm2022; Wu & Lele Reference Wu and Lele2025b )

(3.5)

\begin{align} \tau _{\textit{ij}} &= c_1 \boldsymbol{I} + c_2 \boldsymbol{\bar {S}} + c_3 \boldsymbol{\bar {S}^2} + c_4 \boldsymbol{\bar {\varOmega }^2} + c_5 \big(\boldsymbol{\bar {S}\bar {\varOmega } - \bar {\varOmega }\bar {S}}\big) + c_6 \boldsymbol{\bar {\varOmega }\bar {S}\bar {\varOmega }} + c_7 \big(\boldsymbol{\bar {S}^2\bar {\varOmega } - \bar {\varOmega }\bar {S}^2}\big) \nonumber\\&\quad + c_8 \big(\boldsymbol{\bar {\varOmega }\bar {S}\bar {\varOmega }^2 - \bar {\varOmega }^2\bar {S}\bar {\varOmega }}\big), \end{align}

where $\boldsymbol{\bar {S}}$ and $\boldsymbol{\bar {\varOmega }}$ denote the filtered strain rate or rotation rate tensor. For completeness, the invariants are listed as follows:

(3.6)

\begin{equation} \left\{tr\big(\boldsymbol{\bar {S}}\big), tr\big(\boldsymbol{\bar {S}^2}\big), tr\big(\boldsymbol{\bar {S}^3}\big), tr\big(\boldsymbol{\bar {\varOmega }^2}\big), tr\big(\boldsymbol{\bar {S}\bar {\varOmega }^2}\big), tr\big(\boldsymbol{\bar {S}^2\bar {\varOmega }^2}\big), tr\big(\boldsymbol{\bar {S}^2\bar {\varOmega }^2\bar {S}\bar {\varOmega }}\big) \right\}\!. \end{equation}

The invariants of the tensor basis expansion, which are costly to compute, are substituted with the invariants of the resolved velocity gradient tensor ( $\bar {P}, \bar {Q}, \bar {R}$ ), as well as the magnitudes of the strain rate and rotation rate tensor as inputs instead ( $|\bar {S}|, |\bar {\varOmega }|$ ). The neural network architecture is shown in figure 2, which has many similarities to Wu & Lele (Reference Wu and Lele2025a ).

Figure 2.

Overall GNN architecture. Note that ‘res. conn.’ denotes residual connections, and these are used in conjunction with skip connections to prevent oversmoothing. There are two layer normalisation layers placed in this architecture. Here ReLU is rectified linear unit.

Diagram of a neural network architecture for flow prediction.

As seen, the tensor basis provides the tensorial structure for the data-driven subgrid model, and serves as an inductive bias for the neural network so that any output of the neural network obeys the structure of $\tau _{\textit{ij}}$ . Meanwhile, the inputs (such as $\bar {P}, \bar {Q}, \bar {R}$ ) pass through layers with learnable weights and nonlinear activation functions. These are the inputs that go in from the left-most side of neural network architecture in the diagrams, and will contain features that the neural network will extract and ‘learn’ from. In general, the neural network is more sensitive to inputs that go through nonlinear transformations (Novak et al. Reference Novak, Bahri, Abolafia, Pennington and Sohl-Dickstein2018). Simplifying these inputs (less higher powers of the velocity gradient tensor) would reduce their susceptibility to distribution shifts, leading to more consistent a posteriori performance. To investigate this, neural networks are trained with two different sets of inputs, $\bar {P}, \bar {Q}, \bar {R}, |\bar {S}|, |\bar {\varOmega }|$ (‘complex inputs’) and $\bar {P}, |\bar {S}|, |\bar {\varOmega }|$ (‘simple inputs’). The simpler set is hypothesised to perform more consistently a posteriori. The neural network outputs consist of the coefficients of the tensor basis expansion, which are then multiplied by the tensor basis to obtain $\tau _{\textit{ij}}$ .

The neural network model inputs and outputs should be dimensionless. A global normalisation (max-min normalisation) is used for the input fields

(3.7)

\begin{equation} x_{\textit{norm}} = \frac {x-x_{\textit{min}}}{x_{\textit{max}}-x_{\textit{min}}}*2-1, \end{equation}

while a local normalisation is adopted for the tensor basis and the subgrid stress tensor

(3.8)

\begin{equation} W^{norm}_{\textit{ij}} = \frac {W_{\textit{ij}}}{|\bar {V}_{\textit{ij}}|^x}, \quad \tau _{\textit{ij}}^* = \frac {\tau _{\textit{ij}}}{\varDelta _{\textit{LES}}^2|\bar {V}_{\textit{ij}}|^2}, \end{equation}

where $W_{\textit{ij}}$ is the tensor to be normalised, $\bar {V}_{\textit{ij}}$ denotes the velocity gradient tensor, $|{\bar {V}_{\textit{ij}}}|$ its magnitude, $x$ the exponent to ensure dimensionlessness, $\varDelta _{\textit{LES}}$ corresponds to the LES grid spacing and repeated indices do not imply summation.

To ensure that the results are consistent across slight perturbations to the neural network architecture, a variant architecture is introduced with additional layer normalisation (ALN) operations between select GNN layers, which is shown in figure 3. Input normalisation is also varied, using local normalisation. Training follows Wu & Lele (Reference Wu and Lele2025b ) with a composite loss function that minimises the root mean squared error between the predicted and actual subgrid stress tensor (denoted as non-bolded RMSE when only reporting the value), and the RMSE between the predicted and actual energy dissipation, $\epsilon = - \tau _{\textit{ij}} \bar {S}_{\textit{ij}}$ (denoted as non-bolded DRMSE when only reporting the value):

(3.9)

\begin{equation} \mathcal{L} = \boldsymbol{RMSE}(\tau _{\textit{ij,pred}}-\tau _{\textit{ij}})/\boldsymbol{RMS}(\tau _{\textit{ij}})+ \boldsymbol{RMSE}(\epsilon _{\textit{pred}}-\epsilon )/\boldsymbol{RMS}(\epsilon ) \end{equation}

where $\boldsymbol{RMSE}$ is the root-mean-squared error operation and $\boldsymbol{RMS}$ is the root-mean-square operation.

Figure 3.

Graph neural network variant architecture. Note that ‘res. conn.’ denotes residual connections, and these are used in conjunction with skip connections to prevent oversmoothing. There are five layer normalisation layers placed in this architecture, three more than the original. Other than that, the architecture is exactly the same.

Diagram of a graph neural network variant architecture.

3.2. Training data

3.2.1. Training configuration

The training data are obtained from querying the Johns Hopkins Turbulence Database for forced homogeneous isotropic turbulence (HIT) and channel flow data (Perlman et al. Reference Perlman, Burns, Li and Meneveau2007; Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008). In all cases, the data are split into a 80–10–10 split, where 80 % of the training data are used for training, 10 % of the data are used as a validation set and 10 % of the data are used as the hold-out test set. For HIT, 5 snapshots are queried, and for channel flow, 4 snapshots are queried. Each snapshot is queried at equal time intervals that are statistically de-correlated from each other. The DNS HIT resolution queried is $1024^3$ , with a Taylor Reynolds number ( $Re_{\lambda }$ ) of 433. The channel flow resolution queried is $2048 \times 512 \times 1536$ , with a friction Reynolds number ( $Re_{\tau }$ ) of 1000.

In all cases, the training data are filtered using two different filter widths, 16 $\varDelta _{\textit{DNS}}$ and 32 $\varDelta _{\textit{DNS}}$ , as filtering on two different filter widths shows more generalisability to other filter widths, as seen in Park & Choi (Reference Park and Choi2021). Furthermore, two different training data domain sizes (such as $16^3$ and $32^3$ ) are also used to sensitise the neural network to different domain sizes, following Wu & Lele (Reference Wu and Lele2025b ), where detailed training data processing procedures can be found. This means that the filtered DNS fields will be broken up into ‘subdomains’ where each subdomain is a training data snapshot. If the Kolmogorov scale is defined as $\eta$ , then the total training set for each filter consists of the following, which is consistent with the training set used in Wu & Lele (Reference Wu and Lele2025b ):

(i) 320 realisations of filter width $\bar {\varDelta } = 16\varDelta _{\textit{DNS}}$ ( $\bar {\varDelta }/\eta \approx 35$ ), domain size $16^3$ HIT;
(ii) 320 realisations of filter width $\bar {\varDelta } = 32\varDelta _{\textit{DNS}}$ ( $\bar {\varDelta }/\eta \approx 70$ ), domain size $8^3$ HIT;
(iii) 40 realisations of filter width $\bar {\varDelta } = 16\varDelta _{\textit{DNS}}$ ( $\bar {\varDelta }/\eta \approx 35$ ), domain size $32^3$ HIT;
(iv) 40 realisations of filter width $\bar {\varDelta } = 32\varDelta _{\textit{DNS}}$ ( $\bar {\varDelta }/\eta \approx 70$ ), domain size $16^3$ HIT;
(v) 320 realisations of filter width $\bar {\varDelta } = 16\varDelta _{\textit{DNS}}$ , domain size $16^3$ channel flow;
(vi) 320 realisations of filter width $\bar {\varDelta } = 32\varDelta _{\textit{DNS}}$ , domain size $8^3$ channel flow;
(vii) 40 realisations of filter width $\bar {\varDelta } = 16\varDelta _{\textit{DNS}}$ , domain size $32^3$ channel flow;
(viii) 40 realisations of filter width $\bar {\varDelta } = 32\varDelta _{\textit{DNS}}$ , domain size $16^3$ channel flow.

For each filter, one can generate the above set of training data using the filtered DNS fields. The final training set consists of combining all the entries in the above list for each filter into a large training set. To demonstrate applicability to channel flow, select neural networks are applied to both HIT and channel flow without retraining. The same applies to when two different grid resolutions are used, as some neural networks are applied to a coarser grid spacing to test applicability of the results to a different grid resolution. Even if the neural network trained is only applied to HIT at one specific resolution (such as $\varDelta _{\textit{LES}} = 16 \varDelta _{\textit{DNS}}$ ), they are still trained on the full set of training data with various domain sizes, filter widths and two different flows for consistency.

3.2.2. Effect of filtering on ground truth

Before training the neural networks, an analysis of the ground truth used for supervised learning is presented. Here, the impact of filtering is investigated for the ground truth labels, where both spatial distributions and probability density functions (PDFs) are shown. For the spatial distributions shown, a random slice out of the three-dimensional field was selected for visualisation, while the PDF is plotted for the full three-dimensional field. The ground truth values are shown for the box, BTF and DSCFs.

As seen in figures 4(a) and 4(b) , the act of applying two-thirds dealiasing completely changes the subgrid stress field, which is seen visually in the increasing intensity of the subgrid stress and the additional oscillations in the overall field. Meanwhile, the box and DSCF fields look more similar. This is further seen in the PDF, where the act of two-thirds dealiasing creates far more variance in the PDF as compared with the box filter or the DSCF. The phenomenon happens for both the normal and the shear components of the subgrid stress tensor. This suggests that, when considering the subgrid stress, a Fourier-spectral code with two-thirds dealiasing applied will likely produce a completely different environment than a sixth-order compact code, and the differences can be quite large. Therefore, these two numerical methods are sufficiently different when considering the subgrid stress.

Figure 4.

Effect of various filters on the ground truth subgrid stress quantities. All colour bars are scaled the same per row: (a) $\tau _{\textit{11}}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations; (b) $\tau _{12}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations; (c) $\epsilon = -\tau _{\textit{ij}}\bar {S}_{\textit{ij}}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations.

Multiple graphs depict the spatial distribution and probability density function of subgrid stress quantities under different filters.

Figure 4. Long description

Panel A: Three heatmaps and one line graph. The heatmaps show the spatial distribution of τ11 under Box, BTF, and DSCF filters. The x and y axes correspond to coordinates in the physical domain. The color scale ranges from −0.04 to 0.04. The line graph shows the probability density function of τ11 for the three filters, with the x-axis labeled τ11 and the y-axis labeled Probability density. Panel B: Three heatmaps and one line graph. The heatmaps show the spatial distribution of τ12 under Box, BTF, and DSCF filters. The x and y axes correspond to coordinates in the physical domain. The color scale ranges from −0.04 to 0.04. The line graph shows the probability density function of τ12 for the three filters, with the x-axis labeled τ12 and the y-axis labeled Probability density. Panel C: Three heatmaps and one line graph. The heatmaps show the spatial distribution of ε = −τijSij under Box, BTF, and DSCF filters. The x and y axes correspond to coordinates in the physical domain. The color scale ranges from −0.1 to 0.1. The line graph shows the probability density function of ε = −τijSij for the three filters, with the x-axis labeled −τijSij and the y-axis labeled Probability density.

The same effects are seen in figure 4(c) in the interscale energy transfer field ( $\epsilon = -\tau _{\textit{ij}}\bar {S}_{\textit{ij}}$ ), where the box filter and DSCF produce more similar fields, while the BTF produces far more oscillatory fields. This is reflected in the PDF, where the BTF once again produces far more variance in the PDF as compared with the other filters. Overall, for all the PDFs, the box filter has the largest peak, followed by the DSCF, and then there is quite a large gap between those two filters and the BTF. This also reveals that if one were to train on only one filter (oftentimes a box filter), if applied to a Fourier-spectral code with two-thirds dealiasing, there will automatically be a distribution shift as the neural network is trained to produce subgrid stress distributions similar to the training data, which are filtered with a box filter. As such, training on two or more filters allows the neural network to learn a wider range of plausible subgrid stress distributions it may see in the flow solver, reducing the distribution shift.

3.3. All neural network configurations

A summary of all the neural networks trained are seen in table 1. Neural networks are trained with different training data (one-filter versus augmented training data), different inputs (complex versus simple) and slightly different architectures (more or less layer normalisation). The training spans a wide range of filters, from purely training on only the box filter, only the BTF, only the DSCF, to combinations of two filters and also three filters. Neural networks trained on only the box filter are applied to both the Fourier-spectral with two-thirds dealiasing solver and the sixth-order compact solver. The neural network trained only on the BTF filter is applied to the Fourier-spectral solver to give a baseline on how a neural network trained with only BTF would do. Similarly, the neural network trained with only the DSCF is applied to the sixth-order compact solver to give a baseline. Meanwhile, the neural networks trained with two filters are applied to different flow solvers depending on which combination of filters they were trained on. For example, a neural network trained on the box filter and BTF is applied to the Fourier-spectral code with two-thirds dealiasing (and not to the sixth-order compact code which does not have two-thirds dealiasing). Similarly, a neural network trained on the box filter and DSCF is applied to the sixth-order compact code. The neural network trained with three filters (box, BTF and DSCF) is applied to both solvers without retraining to show how the data augmentation approach can work on different flow solvers without retraining. In all cases, the neural networks used have roughly 24 000 trainable parameters.

Table 1.

Neural network configurations; ‘S’ stands for the spectral solver with two-thirds dealiasing (PadeOps), and ‘C’ stands for the sixth-order compact solver (PadeLibs). Normalisation and neural network structure are also given.

A table comparing neural network configurations with different filters, inputs, and solvers.

Table 1. Long description

The table presents a comparison of neural network configurations. It includes columns for the number of filters, inputs, neural network structure, flow solver, normalization, and neural network name. The table has 18 rows and 6 columns. Row 1: One filter (box), Original, S, C, Global, NN-Box-Complex-Original-G. Row 2: One filter (BTF), Original, S, Global, NN-BTF-Complex-Original-G. Row 3: One filter (DSCF), Original, C, Global, NN-DSCF-Complex-Original-G. Row 4: Two filter (box and BTF), Original, S, Global, NN-BoxBTF-Complex-Original-G. Row 5: Two filter (box and DSCF), Original, C, Global, NN-BoxDSCF-Complex-Original-G. Row 6: Three filter (box, BTF, DSCF), Original, S, C, Global, NN-BoxBTFDSCF-Complex-Original-G. Row 7: One filter (box), Additional LN, S, C, Global, NN-Box-Complex-ALN-G. Row 8: Two filter (box and BTF), Additional LN, S, Global, NN-BoxBTF-Complex-ALN-G. Row 9: Two filter (box and DSCF), Additional LN, C, Global, NN-BoxDSCF-Complex-ALN-G. Row 10: One filter (box), Original, S, C, Global, NN-Box-Simple-Original-G. Row 11: Two filter (box and BTF), Original, S, Global, NN-BoxBTF-Simple-Original-G. Row 12: Two filter (box and DSCF), Original, C, Global, NN-BoxDSCF-Simple-Original-G. Row 13: One filter (box), Additional LN, S, C, Global, NN-Box-Simple-ALN-G. Row 14: Two filter (box and BTF), Additional LN, S, Global, NN-BoxBTF-Simple-ALN-G. Row 15: Two filter (box and DSCF), Additional LN, C, Global, NN-BoxDSCF-Simple-ALN-G. Row 16: One filter (box), Original, S, C, Local, NN-Box-Simple-Original-L. Row 17: Two filter (box and BTF), Additional LN, S, Local, NN-BoxBTF-Simple-ALN-L. Row 18: Two filter (box and DSCF), Additional LN, C, Local, NN-BoxDSCF-Simple-ALN-L.

4. A priori analysis

4.1. One-filter versus two-filter results

As seen in table 2, where the neural networks are evaluated on the test set box filtered data, there is no a priori performance degradation when training on two or more filters versus training on one filter since all RMSE, DRMSE and correlation coefficient values are within 5 % of each other. For example, neural networks trained with a max-min normalisation and using the set of inputs $\bar {P}, \bar {Q}, \bar {R}, |\bar {S}|, |\bar {\varOmega }|$ all have RMSE values in the 0.75–0.78 range, DRMSE values in the 0.66–0.70 range and correlation coefficients in the 0.63–0.66 range. Also, using a complex or simple set of inputs has negligible difference on the neural network performance, suggesting that just using $|\bar {S}|, |\bar {\varOmega }|$ is sufficient ( $\bar {P} = 0$ for incompressible flow, and is manually masked to be zero to prevent the neural network from overfitting to noise. It is kept so that the model can generalise to compressibility in future work).

Table 2.

A priori evaluation on the hold-out test set, each value is given as mean (standard deviation). For the RMSE and DRMSE, lower values are better, while higher values are better for the correlation.

A table comparing the performance of different neural networks on a test set using RMSE, DRMSE, and correlation metrics.

Table 2. Long description

The table presents an a priori evaluation of various neural networks on a hold-out test set, with each value given as mean (standard deviation). The table has 11 rows and 4 columns. The columns are labeled ‘Neural network‘, ‘RMSE’, ‘DRMSE’, and ‘Correlation’. The rows list different neural network configurations and their corresponding performance metrics. Notable trends include similar performance across different neural networks, with RMSE values ranging from 0.750 to 0.772, DRMSE values ranging from 0.661 to 0.693, and correlation coefficients ranging from 0.636 to 0.661. The data suggests that training on two or more filters versus training on one filter does not result in performance degradation, and using a complex or simple set of inputs has negligible difference on the neural network performance.

5. A posteriori analysis

The neural networks are evaluated on two different flow solvers running forced HIT on a $64 \times 64 \times 64$ grid (corresponding to 16 $\varDelta _{\textit{DNS}}$ ) and at $Re_{\lambda } = 433$ . A quantitative metric is also given, sum spectral error ( $SSE = \sum _{k} |(log(E_{DNS_k})-log(E_{LES_k}))|, k_{\textit{PadeLibs}} \in [0, 35], k_{\textit{PadeOps}} \in [0, 24]$ ), where $E_{DNS_k}$ denotes the filtered DNS spectra, $E(k)$ denotes the LES energy spectra, with differences in ‘k’ since PadeOps has two-thirds dealiasing. While SSE is a good first-order metric, it must be interpreted in a physical context. Here, SSE is computed using the logarithm of the energy spectrum values, rendering it artificially sensitive to deviations in the higher-wavenumber, dissipative range. This can manifest in SSE oftentimes assigning a larger SSE value to dissipative spectra even if it is physically desirable, while not putting enough emphasis on high-wavenumber energy buildup. Therefore, visual inspection of the spectral decay or pile-up is still essential. Other metrics that are investigated for HIT include plotting the PDF of the production term of the enstrophy equation, $\bar {\omega }_i \bar {S}_{\textit{ij}} \bar {\omega }_{\!j}$ , as well as analysing the intermittency by plotting longitudinal flatness as a function of spatial spacing, where longitudinal flatness is defined as

(5.1)

\begin{equation} F(r) = \frac {\langle (\delta u)^4 \rangle }{\langle (\delta u)^2 \rangle ^2} \end{equation}

and $\delta u$ is defined as $\delta u = \bar {u}(x+r) - \bar {u}(x)$ , where r is the spatial spacing.

5.1. Fourier-spectral with two-thirds dealiasing (PadeOps)

The neural networks are integrated into the Fourier-spectral code with two-thirds dealiasing (PadeOps), and after running to statistical stationarity, data are collected and plotted. In addition, two traditional models are run for reference. The static Smagorinsky model is used with a coefficient set to $C = 0.18$ , while the Sigma model is also run with a coefficient set to $C=1.3$ . These are all within the ranges commonly used for these two models (Smagorinsky Reference Smagorinsky1963; Meyers & Sagaut Reference Meyers and Sagaut2006; Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011). The Sigma model was chosen as the coefficient range works for many different flows and does not require clipping for channel flow, whereas models such as the dynamic Smagorinsky model may need clipping (Meneveau et al. Reference Meneveau, Lund and Cabot1996). Table 3 shows all the neural networks trained and integrated into the Fourier-spectral solver, as well as the SSE for easy reference.

Table 3.

A posteriori evaluation for the Fourier-spectral code for all cases tested.

A table comparing the performance of various neural networks and traditional models based on their sum of squared errors (SSE).

Table 3. Long description

The table presents a comparison of the sum of squared errors (SSE) for various neural networks and traditional models. It consists of two columns and multiple rows. The first column lists different neural networks, and the second column shows their corresponding SSE values. The neural networks include NN-Box-Complex-Original-G, NN-BTF-Complex-Original-G, NN-BoxBTF-Complex-Original-G, NN-BoxDSCFBTF-Comp-Orig-G, NN-Box-Complex-ALN-G, NN-BoxBTF-Complex-ALN-G, NN-Box-Simple-Original-G, NN-BoxBTF-Simple-Original-G, NN-Box-Simple-ALN-G, NN-BoxBTF-Simple-ALN-G, NN-Box-Simple-Original-L, and NN-BoxBTF-Simple-Original-L. The traditional models listed are Static Smagorinsky and Sigma. The SSE values for the neural networks range from 3.68 to 18.87, while the SSE values for the traditional models are 16.6 for Static Smagorinsky and 13.34 for Sigma. The table provides a clear comparison of the performance of different neural networks and traditional models based on their SSE values.

5.1.1. Effect of multiple filters

First, the effect of training on two filters as compared with one filter is analysed, followed by the analysis of a three-filter approach using the ‘complex’ set of inputs. Figures 5(a) to 5(c) show the results of the neural network trained only with a box filter and the neural network trained with two filters (box and BTF) compared with the static Smagorinsky and Sigma models.

Figure 5.

A posteriori one- and two-filter comparisons with traditional models for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Three line graphs compare different models’ performance in Fourier-spectral code with two-thirds dealiasing.

From the spectra, it is seen that the static Smagorinsky and Sigma models are over-dissipative, despite using the smallest coefficient in the recommended range for the Sigma model. This highlights that in the coarse-grid regime, traditional models may be slightly less accurate as large amounts of turbulent kinetic energy are unable to be resolved and the subgrid stress models must model larger fractions of the unresolved physics. Meanwhile, the neural network trained on one filter is under-dissipative, causing an unrealistic buildup of turbulent kinetic energy at the higher wavenumbers. The neural network trained on two filters is most able to accurately represent the spectra. The spectral fall-off for all of the LES results is expected when compared with filtered DNS due to the two-thirds dealiasing operation at the higher wavenumbers. When analysing the longitudinal flatness, which characterises the intermittency at different spatial scales, one can see that the two-filter case is more accurate than the one-filter case. However, in both these cases the neural network models underpredict the filtered DNS values, while the traditional models are more accurate in representing these values. This is expected, as over-dissipative models tend to reduce the variance (denominator of the longitudinal flatness calculation), resulting in higher values of ‘rare events’ compared with the variance. However, an intermittent flow that does not drive turbulence for HIT would not be as physically accurate.

One of the primary drivers of three-dimensional forced HIT is vortex stretching, which can be analysed by plotting the production term of the resolved enstrophy equation. From figure 5(c), the neural network models are revealed to drive enstrophy production far more accurately than the traditional models. The traditional models ‘fall off’ faster than the filtered DNS values, while the neural network models are more accurate, suggesting that they are driving the energy cascade in a more physical manner. The two-filter case is also more accurate than the one-filter case, which can be seen in the positive tail of the enstrophy production. The one-filter case, which over-predicts the enstrophy production, when combined with the low longitudinal flatness values at small separation scales, suggests that there is grid-scale noise that is not being adequately dissipated, which is confirmed in the high-wavenumber buildup in the spectra. Overall, when combining all this analysis, the two-filter case is far more physically accurate than the one-filter case, and in general the two-filter neural network case is more physically accurate than the traditional models in representing turbulent physics. This suggests that filtering plays a crucial role in the accuracy of neural network subgrid stress models that are trained outside of the LES solver. Oftentimes, for the enstrophy production, the PDF tails are lower when compared with filtered DNS (less variance). This can partially be attributed to the spectral fall-off induced by the two-thirds dealiasing scheme, as enstrophy production is dominated by the highest resolved wavenumbers. Thus, the two-thirds dealiasing scheme will cause slightly smoothed velocity gradients, which will lead to smaller peak magnitudes of enstrophy production and a tail that ‘falls off’ faster in the PDF.

Next, the box-filter-trained neural network and the two-filter neural network are compared with a neural network trained on all three filters (box, BTF, DSCF) and a neural network trained on only the BTF filter. This will highlight the neural network performance when trained only on BTF, as well as the neural networks trained on all three filters, which can be applied without retraining to the sixth-order compact LES solver too. The total comparison is shown in figures 6(a) to 6(c).

Figure 6.

A posteriori multiple filter comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation, but trained with a variety of filter combinations ranging from one to three different filter shapes.

Multiple graphs compare neural network performance using different filter combinations.

Figure 6. Long description

Panel A: A log-log line graph displays energy spectra E(k) against wavenumber k. The graph includes multiple lines representing different neural network models and a filtered DNS. The legend indicates the models and their respective sum of squared errors (SSE). Panel B: A line graph shows longitudinal flatness F(r) as a function of separation scale r. Different lines represent various neural network models and filtered DNS, with the legend specifying each model. Panel C: A log-scale probability density function (PDF) graph depicts the distribution of a variable against its values. Multiple lines represent different neural network models and filtered DNS, with the legend identifying each model.

The spectra indicate that the neural network trained with only the BTF filter is more accurate than the neural network trained with only the box filter, which is to be expected as BTF accounts for two-thirds dealiasing while the box filter does not. Meanwhile, the three-filter neural network also does well in producing flows with physically accurate spectra. As seen in the SSE from table 3 and figure 6(a), the SSE is the largest for the box-filter case, while the other neural networks trained have more reasonable SSE. However, when taking a look at the longitudinal flatness, the benefits of training on multiple filters are still clear. For the neural network trained on just the BTF, at small separation scales, the overall intermittency of the flow is still lower as compared with both the two- and three-filter cases. In fact, the neural networks trained on just one filter (whether it is box or BTF) are the worst, while the multiple-filter cases are better in capturing the intermittency. Overall, all cases have similar behaviour at larger separation scales, although the three-filter case overpredicts the longitudinal flatness, but is also the best at predicting the flow intermittency at the smallest separation scales. When analysing the production of enstrophy, the BTF-trained neural network and the three-filter-trained neural network show similar behaviour, although the three-filter-trained neural network is slightly more accurate (closer to filtered DNS). For this case, the two-filter (box and BTF) trained neural network is still the most accurate. In general, when considering the spectra, intermittency and the production of enstrophy, the two cases that are closest to filtered DNS are the two- and three-filter-trained neural networks. This suggests that it is beneficial for the neural network to be trained on multiple filter shapes, as opposed to one filter shape. By doing so, the neural network is less likely to ‘overfit’ to a specific filter shape and can become a more robust model for a posteriori LES applications, where the implied filter is unknown. The implied filter width for LES is a product of both the grid and numerical method, and thus the implied SGS distribution is entirely dependent on the implied filter. By training on multiple different filter shapes, the output distribution shift for $\tau _{\textit{ij}}$ can be better reduced.

5.1.2. Neural network perturbation

The box-filtered and two-filter training data are also used to train a neural network variant with ALN to evaluate if the accuracy increase due to the two-filter training is present across neural network perturbations. These results are shown in figures 7(a) to 7(c).

Figure 7.

A posteriori ALN comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the variant ALN architecture with the global normalisation.

Three line graphs compare different data sets.

From the spectra, the one-filter-trained neural network is now significantly better as compared with the previous case, suggesting that the results are very sensitive to the neural network architecture. Meanwhile, the two-filter results are more robust across a wider range of neural network architectures, as the shape of the spectra does not change much. As seen, there is still an improvement in the spectra when training on two filters as compared with one filter, although the effect is less drastic for this neural network architecture perturbation. However, the longitudinal flatness still suggests that the two-filter case is better at representing the overall ‘tails’ of the flow, as it shows better small-scale intermittency. Meanwhile, both neural network models are accurate in enstrophy production, with the one-filter case producing a more accurate PDF at the negative tail but the two-filter case producing a more accurate PDF at the positive tail. Overall, the two-filter case is still more physically accurate as compared with the one-filter case, and when comparing with the previous case, the two-filter case is far more robust to architecture perturbations than the one-filter case.

5.1.3. Using simpler inputs

Now, neural networks are trained on the simpler inputs to evaluate the distribution shift associated with neural network inputs. The hypothesis is that simpler inputs suffer less from a distribution shift, and so there should be less variation in these results. The neural networks trained involve both the variant architecture and the original architecture, as well as a different input normalisation to show that the results are robust to a wide range of conditions. The results are shown in figures 8(a) to 8(c).

Figure 8.

A posteriori simple input comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. Both the original and ALN GNN architectures are tested with ‘simple’ inputs using both local and global normalisation.

Three line graphs compare different data sets related to neural network performance in subgrid stress modeling.

Figure 8. Long description

Panel A: A line graph displays the energy spectrum E(k) against the wavenumber k. The graph includes multiple lines representing different neural network models and a filtered DNS. The legend indicates various models with their respective structure similarity error (SSE) values. Panel B: A line graph shows the longitudinal flatness F(r) as a function of the separation scale r. Multiple lines represent different neural network models and a filtered DNS, with a legend detailing each model. Panel C: A line graph presents the probability density function (PDF) on a log scale against the normalized variable. Multiple lines represent different neural network models and a filtered DNS, with a legend specifying each model.

Across all the different neural network models, the spectra are all relatively consistent, which validates the idea that simple inputs result in less distribution shift as compared with more complex inputs. From the SSE from table 3 and figure 8(a), the neural networks trained with simple inputs have relatively similar SSE values across a wide range of training conditions. When evaluating the longitudinal flatness at small separation scales, the two-filter-trained neural networks in general, predict a higher (closer to filtered DNS) value as compared with the neural networks trained on one filter. The enstrophy production PDFs show relatively similar performance across all the various models. This suggests that training on two filters will often do no worse than training on one filter, and can still obtain some performance improvements in select situations, highlighting the utility of training on two filters. This is seen in two different neural network architectures and two different input normalisations. Overall, using simpler inputs does reduce the distribution shift seen a priori versus a posteriori, and training on one versus two filters has less of a difference, but can still improve some qualities of the flow. Ultimately, this suggests that the input distribution shift has a large effect on the neural network SGS models a posteriori, as simpler inputs (which suffer less from distribution shift) produce neural network models that have less a posteriori variance across models. When the input distribution shift is reduced, training on two filters can still help increase robustness a posteriori by reducing the output distribution shift, although this becomes a more secondary effect.

5.2. Sixth-order compact (PadeLibs)

The neural networks are now integrated into a sixth-order compact code, PadeLibs, to investigate if the results hold for a flow solver with a different numerical method. Similarly, forced HIT is run with the same parameters as the Fourier-spectral code, and the spectra, longitudinal flatness and enstrophy production PDFs are plotted. Table 4 shows all the neural networks trained and integrated into this flow solver, as well as SSE for easy reference.

Table 4.

A posteriori evaluation for the sixth-order compact code for all cases tested.

A table comparing neural networks and their SSE values for a sixth-order compact code.

Table 4. Long description

A table with two columns and eleven rows comparing different neural networks and their corresponding SSE values. The first column lists neural networks, and the second column lists their SSE values. Row 1: NN-Box-Complex-Original-G, 32.58. Row 2: NN-DSCF-Complex-Original-G, 35.33. Row 3: NN-BoxDSCF-Complex-Original-G, 11.93. Row 4: NN-BoxDSCFBTF-Comp-Orig-G, 14.22. Row 5: NN-Box-Complex-ALN-G, 9.50. Row 6: NN-BoxDSCF-Complex-ALN-G, 5.59. Row 7: NN-Box-Simple-Original-G, 7.56. Row 8: NN-BoxDSCF-Simple-Original-G, 10.97. Row 9: NN-Box-Simple-ALN-G, 9.50. Row 10: NN-BoxDSCF-Simple-ALN-G, 12.22. Row 11: NN-Box-Simple-Original-L, 12.56. Row 12: NN-BoxDSCF-Simple-Original-L, 15.84.

5.2.1. Effect of multiple filters

First, the neural network trained on only the box filter is analysed against the two-filter case. In this case, the two-filter training data correspond to data prepared with both the DSCF and the box filter. These results are shown in figures 9(a) to 9(c).

Figure 9.

A posteriori one- and two-filter comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Three line graphs compare different data sets across various metrics.

Figure 9. Long description

Panel A: A line graph displays the energy spectrum E(k) against the wavenumber k. The x-axis is labeled ‘k’ and the y-axis is labeled ‘E(k)’, both on a logarithmic scale. Three data sets are shown: Filtered DNS in black, NN-Box-Complex-Original-G in blue, and NN-BoxDSCF-Complex-Original-G in red. The SSE values for the blue and red lines are 32.58 and 11.93, respectively. Panel B: A line graph shows the longitudinal flatness F(r) against the separation scale r. The x-axis is labeled ‘Separation scale r’ and the y-axis is labeled ‘Longitudinal flatness F(r)’. Three data sets are shown: Filtered DNS in black, NN-Box-Complex-Original-G in blue, and NN-BoxDSCF-Complex-Original-G in red. Panel C: A line graph presents the probability density function (PDF) on a logarithmic scale against the variable 〈ωiSijωj〉. The x-axis is labeled ‘〈ωiSijωj〉’ and the y-axis is labeled ‘PDF (log scale)’. Three data sets are shown: Filtered DNS in black, NN-Box-Complex-Original-G in blue, and NN-BoxDSCF-Complex-Original-G in red.

The spectra indicate that the neural network trained on two filters is more physically accurate than the neural network trained on one filter, which suffers from an extreme buildup of high-wavenumber energy. This is consistent with the previous Fourier-spectral results, showing robustness across two different high-order numerical methods. In addition, the small-scale longitudinal flatness also shows improvements when training on two filters instead of one filter, which is also consistent with the previous results on the Fourier-spectral code. For the large-scale longitudinal flatness, both cases flatten out at higher values as compared with the filtered DNS profile. The enstrophy production PDF plot shows that the neural network trained on one filter is significantly over-predicting the enstrophy production and, combined with the smaller longitudinal flatness value at small separation scales (suggests high variance compared with the intermittent events), highlights that there is spurious grid-scale energy that is not being dissipated by the model. Overall, the results indicate that training on two filters is far better than training on only a box filter. This suggests that for a solver with a different high-order numerical method, the idea of training on two filters to reduce the output distribution shift of $\tau _{\textit{ij}}$ holds, and produces more robust models a posteriori.

Next, these results are contrasted with a neural network trained on only the DSCF, as well as the same three-filter neural network that was used in the Fourier-spectral solver. These results are shown in figures 10(a) to 10(c).

Figure 10.

A posteriori multiple-filter comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation, but trained with a variety of filter combinations ranging from one to three different filter shapes.

The neural network trained with just the DSCF is over-dissipative, which highlights the utility of training on two or more filters. In this case, while the three-filter case is slightly more under-dissipative as compared with the two-filter case, it is still far better than either of the one-filter cases. This can also be seen in table 4 and figure 10(a), where the SSE values for the two-filter and three-filter cases are lower than either of the one-filter cases. The DSCF-trained neural network produces a longitudinal flatness curve that is more similar to the dissipative traditional models. This is expected, as the DSCF-trained neural network is over-dissipative too, and in this case actually over-predicts the small-scale longitudinal flatness. When combined with the over-dissipative nature of DSCF, this suggests that the overall ‘background’ variance of the DSCF-trained neural network simulation is lower, which produces a large longitudinal flatness. Meanwhile, the three-filter case in general has a more accurate longitudinal flatness than the single box-filter case and the two-filter case, even at the smaller separation scales. When evaluating the enstrophy production PDFs, the neural networks trained with two or three filters outperform the box-filter-trained neural network. The DSCF-trained neural network produces an enstrophy production PDF that is also reasonable. When compared with the two-filter case, it is better at predicting the negative tail but is worse at predicting the positive tail, although both underpredict the PDF. However, the DSCF achieves this prediction by sacrificing the spectra, as the rounded shape shows little indication of an inertial range, suggesting that it is unable to capture the overall physics. Meanwhile, the three-filter case over-predicts the PDF at both the positive and negative tails, a potential consequence of training on the oscillatory BTF when the flow solver does not apply two-thirds dealiasing. In general, the two- and three-filter cases retain more of an inertial range in the energy spectrum, while also doing well in predicting the enstrophy production PDFs. Overall, this suggests that training on two or three filters results in models that are more physically accurate than training on one filter, and that the three-filter case produces accurate results across both flow solvers, although it may not be the most optimal model as it was trained on a huge variety of filters, not all of which are applicable to the flow solver (such as training on BTF when the flow solver does not use two-thirds dealiasing). Despite this, three-filter training can be applied to solvers with different numerical methods without re-training, and outperforms training on one filter. Overall, this suggests that training on multiple filters to reduce the neural network output distribution shift of $\tau _{\textit{ij}}$ while ensuring that the neural network model does not ‘overfit’ to the specific characteristics of one specific filter shape is an effective way of training more robust models a posteriori.

5.2.2. Neural network perturbation

Similar to the Fourier-spectral code, the box-filter training data and the two-filter training data (box filter and DSCF) are used to train neural networks with the variant architecture. These results are shown in figures 11(a) to 11(c).

Figure 11.

A posteriori ALN comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the ALN architecture with the global normalisation.

Figure 11. Long description

Three line graphs compare different data sets related to a sixth-order compact code. Panel A: A line graph plots energy spectra E(k) against wavenumber k on a log-log scale. The graph includes three lines representing filtered DNS, NN-Box-Complex-ALN-G, and NN-BoxDSCF-Complex-ALN-G. Panel B: A line graph shows longitudinal flatness F(r) against separation scale r. The graph includes three lines representing filtered DNS, NN-Box-Complex-ALN-G, and NN-BoxDSCF-Complex-ALN-G. Panel C: A line graph displays the probability density function (PDF) of a variable on a log scale against the variable itself. The graph includes three lines representing filtered DNS, NN-Box-Complex-ALN-G, and NN-BoxDSCF-Complex-ALN-G.

When evaluating the spectra, the same trend from the Fourier-spectral code also holds, where the neural network trained with one filter has a wide range of results depending on the neural network architecture, while the neural network trained with two filters has less variance in the results depending on neural network architectures. In this case, both neural networks trained predict the spectra accurately. The longitudinal flatness, especially at smaller separation scales, also shows less variance across the various neural network architectures when trained on multiple filters as compared with one filter. This is seen as the longitudinal flatness at small separation scales for the neural network trained on only the box filter varies significantly depending on whether or not the ALN architecture is used (comparing figures 10 b and 11 b). Meanwhile, the enstrophy production PDF shows little difference between the model trained with one or two filters. This highlights that, once again, the general trend of training on two filters is that it often does no worse than training on one filter, and oftentimes increases the robustness of neural network models a posteriori. This trend holds across two different neural network architectures.

5.2.3. Using simpler inputs

In addition, to ensure that simpler inputs do result in a smaller distribution shift for a different numerical method, neural networks are trained and evaluated with the simple set of inputs, across two different neural network architectures and also two different input normalisations. These results are shown in figures 12(a) to 12(c).

Figure 12.

A posteriori simple input comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. Both the original and ALN GNN architectures are used with the global and local normalisation for the ‘simple’ inputs.

When evaluating the spectra, the simpler inputs reduce the input distribution shift, as the neural networks trained only on the box filter are also able to produce spectra that are reasonable. There is little difference in the spectra when comparing them with the filtered DNS curve. This is also seen in the relatively smaller variance of the SSE that is shown in table 4 and also figure 12(a), especially when compared with the larger variance seen in the complex input cases when training on one versus two filters. A similar trend holds for the longitudinal flatness, where the models produce relatively similar results, regardless of training on one or two filters. One exception is the local normalisation, where training on two filters produces better performance than training on one filter. The enstrophy production plot is also relatively similar, with all models doing well. The exception is the neural networks trained with the local normalisation, where they slightly under-predict the tails more than the other models, although for the negative tail, training on two filters does help reduce the behaviour slightly. Overall, using simpler inputs does reduce the distribution shift in a posteriori simulations, and two-filter training with simpler inputs does relatively similarly to the one-filter training. This shows that for simpler inputs, there is no significant performance degradation when training on two filters as compared with one filter.

6. Additional a posteriori evaluations

To ensure that these trends hold for a different grid resolution and for a different flow, a coarser forced HIT is run on a $32 \times 32 \times 32$ grid (corresponding to a filter width of $\varDelta _{\textit{LES}} = 32 \varDelta _{\textit{DNS}}$ ), which is twice as coarse as the $64 \times 64 \times 64$ grid used in the previous HIT simulations. Furthermore, a channel flow simulation is run, matching the $Re_{\tau } = 1000$ of the DNS simulation used for training. The grid used is $128 \times 96 \times 32$ , which corresponds to a filter width of $\varDelta _{\textit{LES}} = 16 \varDelta _{\textit{DNS}}$ . The wall model used is not a data-driven wall model. Instead, a classic wall model is used that is identical to the one described by Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005). A RK3 time stepping scheme is employed, as well as Fourier collocation in the wall-parallel directions and a sixth-order staggered scheme in the wall-normal direction. In all cases, the Fourier-spectral code with two-thirds dealiasing is run as the primary code.

6.1. Evaluation on a different grid width

The first main trend seen was that training on two filters, in general, is better than training on one filter (especially the box filter), and this trend is more prominently seen for complex inputs. The second trend seen was that training on a simple set of inputs results in less of an input distribution shift, resulting in more consistent results across the various neural network models. In addition, the two-filter cases do no worse than the one-filter cases in many situations. These trends are now evaluated on a forced HIT simulation that is twice as coarse as the previous simulation. To illustrate these trends, a minimal amount of simulations were run. For the complex inputs, the box-filter-trained neural network was run, as well as the box- and BTF-trained neural network. This is shown to highlight the performance improvement that can be gained over using a filter that is traditionally used in literature (box filter). These results are shown in figures 13(a) to 13(c). To illustrate the simple input set, a similar filter configuration is run, with results shown in figures 14(a) to 14(c).

Figure 13.

A posteriori one- versus two-filter comparison with the Sigma model as a reference for a coarser grid spacing corresponding to $\varDelta _{\textit{LES}} = 32 \varDelta _{\textit{DNS}}$ . Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Three line graphs compare different models and their performance metrics.

Figure 13. Long description

Panel A: A line graph displays energy spectra E(k) against wavenumber k. The graph includes four lines representing different models: NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, Sigma, and Filtered DNS. The SSE values for NN-Box-Complex-Original-G and NN-BoxBTF-Complex-Original-G are 3.33 and 0.87, respectively. Panel B: A line graph shows longitudinal flatness F(r) against separation scale r. The graph includes four lines representing Filtered DNS, NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, and Sigma. Panel C: A line graph displays the probability density function (PDF) on a log scale against the variable 〈ωiSijωj〉. The graph includes four lines representing Filtered DNS, NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, and Sigma.

Figure 14.

A posteriori simple input comparisons with the Sigma model as a reference on the coarser grids. Spectra, longitudinal flatness and the PDF of $\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation applied to the ‘simple’ inputs.

Three line graphs compare different models on coarser grids.

From figure 13(a), when comparing the spectra of the one-filter versus two-filter case, the two-filter case is able to improve the spectra by increasing the dissipation as compared with the one-filter neural network model, which is under-dissipative. For reference, the Sigma model is also shown, and it is over-dissipative, although less so for this grid resolution, despite using the same coefficient as the previous simulations. The same trends hold for the longitudinal flatness, where the traditional models capture this better, while the two-filter-trained neural network significantly improves upon the one neural network case. For the enstrophy production, the neural network models do relatively similarly, but are able to outperform the traditional models such as the Sigma model. However, all models are unable to accurately capture the PDF of the enstrophy production term, likely due to the very coarse grid resolution, as in the finer resolution simulations, the two-filter neural network model was able to more accurately capture the enstrophy production PDF when compared with the one-filter case. Overall, for the complex set of inputs, especially when comparing the energy spectra and longitudinal flatness, training with two filters is still better than training with one filter, showing that this trend still holds even at a different resolution.

Meanwhile, when considering the simple inputs (figures 14 a to 14 c), the neural networks do not show as much difference in the spectra when training with one or two filters, which is consistent with the previous cases. For longitudinal flatness, there is still improvement when training on two filters instead of one filter, especially at the smaller separation scales. The enstrophy production PDFs look similar for both neural network models, although they are still better than the traditional models. Overall, the same trend here holds, where for the simple inputs, there is less of a difference when training on one versus two filters, and that the two-filter case often does no worse than the one-filter case, and can sometimes do better in some metrics such as longitudinal flatness. This suggests that the mechanisms of reducing input and output distribution shift are not specific to one grid resolution, and are general across a range of grid resolutions. Specifically, using simpler inputs to the neural network can reduce the input-side distribution shift, which the neural network is sensitive to since these are the features the neural network is ‘learning’ from. Meanwhile, training on multiple filters can reduce the output-side distribution shift and ensure that the neural network does not ‘overfit’ to one specific filter shape.

6.2. Evaluation on channel flow

Now, these main trends are also illustrated for channel flow. Once again, the box filter is chosen as the filter used for single-filter training to highlight the performance improvement that can be gained using a filter that is traditionally used in the literature, and a minimal set of simulations were run to illustrate these trends. However, for channel flow, an additional set of neural networks are trained with the variant neural network architecture (ALN) to highlight that the variant neural network architecture, while being more robust, can also be enhanced with the two-filter training as compared with the one-filter training.

6.2.1. Two-filter complex input analysis

First, the box-filter-trained neural network is compared with the box- and BTF-trained neural network, along with the traditional models. The spectra and mean profiles are shown in figures 15(a) to 15(b). The spectra reported is for the streamwise velocity, at a particular wall-normal location of $z/H = 0.46875$ .

Figure 15.

A posteriori streamwise velocity spectra and mean profiles for channel flow, comparing one- versus two-filter GNN models with traditional models. The GNN architecture used corresponds to the original architecture with the global normalisation.

Two line graphs compare different models for channel flow, showing streamwise velocity spectra and mean profiles.

Figure 15. Long description

The image contains two line graphs labeled as Panel (a) and Panel (b). Panel (a) is a line graph depicting the streamwise velocity spectra (E(k)) against the wavenumber (k). The vertical axis represents E(k) on a logarithmic scale ranging from 10^−6 to 10^−3, while the horizontal axis represents the wavenumber (k) on a logarithmic scale ranging from 10^0 to 10^1. The graph includes multiple lines representing different models: NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, Filtered DNS, Static, and Sigma. Each model is distinguished by different colors and line styles, with a legend indicating their respective SSE values. Panel (b) is a line graph showing the mean profiles (U) against the normalized height (z/H). The vertical axis represents the mean velocity (U) ranging from 0.2 to 1.2, while the horizontal axis represents the normalized height (z/H) on a logarithmic scale ranging from 10^-2 to 10^0. Similar to Panel (a), this graph also includes multiple lines representing different models: NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, Filtered DNS, Static Smagorinsky, and Sigma. The lines are distinguished by different colors and styles, with a legend indicating the models used.

In this case, the one-filter neural network completely fails, creating an unphysical spectra that causes the channel flow velocity fields to be extremely choppy. In addition, the static Smagorinsky model also has a large buildup of turbulent kinetic energy at the high wavenumbers. Meanwhile, the Sigma model and the neural network model trained on two filters are able to accurately represent the spectra. This is also seen in the mean profiles, where the single filter neural network has poor agreement with the DNS profile. Meanwhile, the two-filter case lines up very well with the filtered DNS profile, even when compared with the traditional models. This provides early indications that training on two filters is also more robust than training on one filter for channel flow. The Reynolds stress profiles are shown in figures 16(a) to 16(c).

Figure 16.

A posteriori one- versus two-filter Reynolds stress comparisons with traditional models for channel flow. Here, $z/H = 0$ corresponds to the wall while $z/H=1$ corresponds to the channel centreline. The GNN architecture used corresponds to the original architecture with the global normalisation.

Three line graphs compare Reynolds stress components in channel flow.

The Reynolds stress profiles show a similar story, where the one-filter case oftentimes performs the worst, overpredicting the Reynolds stress profiles, especially for $R_{\textit{11}}$ and $R_{33}$ . In addition, the static Smagorinsky model also predicts poor Reynolds stress profiles. Meanwhile, the neural network trained on both the box filter and BTF is able to produce far more accurate Reynolds stress profiles, oftentimes producing the most accurate profiles even when compared with the Sigma model. Improvements due to the one- versus two-filter training is clearly seen in all the Reynolds stress profiles, where the one-filter over-predictions are significantly reduced when training on two filters. Thus, for the complex set of inputs, the overall trend of training on two different filters being more accurate that was seen in forced HIT is also seen in channel flow. This suggests that the idea of reducing the output distribution shift and ensuring that the neural network sees multiple filter shapes through two-filter training is more general as it shows consistent improvement across both forced HIT and channel flow for complex inputs.

6.2.2. Two-filter complex input analysis ALN

Now, the variant architecture (ALN) is analysed and compared with the original architecture for channel flow. The streamwise velocity spectra and streamwise velocity mean profiles are shown in figures 17(a) and 17(b).

Figure 17.

A posteriori spectra and mean profile comparison of different neural network perturbations. Both the original and ALN GNN architectures are used, with global normalisation applied to the ‘complex’ inputs.

Figure 17. Long description

Panel A: A line graph displays the energy spectrum E(k) against the wavenumber k. The x-axis is labeled k and ranges from approximately 1 to 20. The y-axis is labeled E(k) and ranges from 10^−6 to 10^−2. Multiple lines represent different neural network perturbations: NN-Box-Complex-Original-G, NN-BoxBTF-Complex-Original-G, NN-Box-Complex-ALN-G, NN-BoxBTF-Complex-ALN-G, and Filtered DNS. Each line shows different trends in the energy spectrum across the wavenumber range. Panel B: A line graph shows the mean velocity profile 〈U〉 against the normalized height z/H. The x-axis is labeled z/H and ranges from 10^−2 to 1. The y-axis is labeled 〈U〉 and ranges from 0.2 to 1.2. Similar to Panel A, multiple lines represent different neural network perturbations and Filtered DNS, illustrating how the mean velocity profile varies with height.

Overall, both the one-filter and two-filter cases with the variant neural network architecture do well, especially when compared with the one-filter case with the original architecture. This is also shown in the mean profiles, where the majority of the neural network models do well while the original neural network architecture trained with the box-filter data does poorly. From both the spectra and the mean profiles, the training on two filters significantly reduces the variance in the results across neural network architectures as compared with training on one filter. This suggests that, overall, the trend seen from forced HIT holds, where the neural network that is trained with one filter is far more sensitive to neural network architectures, while the two-filter-trained neural networks are more robust across neural network architectures. The Reynolds stress profiles are seen in figures 18(a) to 18(c).

Figure 18.

A posteriori Reynolds stress comparison of different neural network perturbations. Here, $z/H = 0$ corresponds to the wall while $z/H=1$ corresponds to the channel centreline. Both the original and ALN GNN architectures are used, with global normalisation applied to the ‘complex’ inputs.

Three line graphs compare Reynolds stress data across different neural network perturbations.

For $R_{\textit{11}}$ , there is not too much difference for the neural network models that do well. However, for $R_{13}$ and $R_{33}$ , one can see that training on two filters still helps improve the model as compared with training on one filter, even for the variant neural network architecture. In this case, the red curves are more close to the filtered DNS profile as compared with the blue curves. This leads to a more nuanced conclusion for the ALN (variant) neural network architecture.

First, the ALN acts as a stabiliser that prevents the single-filter model from producing wildly inaccurate results, essentially masking its lack of robust physical knowledge that can generalise a posteriori. In addition, the two-filter training provides robustness, since the network trained on two filters learns a more general representation of the physics that is accurate regardless of whether the stabilising ALN layers are present. This means that the two-filter case is robust to architectural perturbations. Even with the ALN neural network architecture, as seen in the Reynolds stress profiles, one can still produce a more physically accurate model when training on two filters versus one. This suggests that by incorporating two filters into the training process, the neural network is able to learn a more robust representation of the SGS physics to produce more accurate turbulence in the LES flow fields, regardless of architecture choice. The differences between one- and two-filter training may be large or small depending on the architecture choice, but there is often improvement in either architecture when training with two filters.

6.2.3. Simple inputs

Now, neural networks are trained with simple inputs and tested in channel flow to evaluate if there is also less input distribution shift for a more complex flow when using simpler inputs. The spectra and mean velocity profiles are shown in figures 19(a) and 19(b).

Figure 19.

A posteriori simple inputs streamwise velocity spectra and mean profiles. The GNN original architecture is used, with global normalisation applied to the ‘simple’ inputs.

Two line graphs compare different data sets related to streamwise velocity spectra and mean profiles.

Overall, similar trends are seen for the spectra and the mean profiles that were observed in forced HIT. In these cases, there is less of a difference between the neural networks trained with one filter or two filters as there is less input distribution shift, so neural networks behave more similarly a posteriori. A similar trend is seen in figures 20(a) to 20(c) for the Reynolds stress profiles. Despite less difference, it can be seen that in some situations, particularly $R_{33}$ , training on two filters produces more physically accurate profiles as compared with training on one filter. This establishes the trend that was seen in forced HIT, that for simpler inputs, there is less of a difference between training with one or two filters, and that the two-filter case often does no worse than the one-filter case, and in some situations can do better.

Figure 20.

A posteriori simple inputs Reynolds stress profiles for channel flow. Here, $z/H = 0$ corresponds to the wall while $z/H=1$ corresponds to the channel centreline. The GNN original architecture is used, with global normalisation applied to the ‘simple’ inputs.

Three line graphs compare Reynolds stress profiles for channel flow.

7. Combined a posteriori analysis

Overall, combining data augmentation and simpler inputs yield robust a posteriori performance across solvers with different high-order numerical methods. Training with two filters allows the neural network to see more diverse input and output distributions, leading to reduced distribution shift a posteriori. This allows the neural network, even though it is trained on generic filters not exactly mimicking the LES implicit transfer function, to perform more accurately to unseen implicit filters as it is not overfitting to the kernel artefacts associated with one specific filter (e.g. Gibbs oscillations). This is seen as the neural network trained on all three filters (and is used in both flow solvers) performs better than the neural network trained only on the box filter (that is also applied to both flow solvers). Meanwhile, the ‘simple’ training set produces results that are more consistent than the more complex training set because it is able to reduce the input distribution shift, as higher-order powers of the velocity gradient tensor suffer more from numerical artefacts such as aliasing, especially at the grid cutoff resolution. By relying on the lower-order terms, the ‘simple’ inputs allow the neural network to be less exposed to and not overfit to high-wavenumber numerical artefacts, increasing the a posteriori robustness of the neural network model. As seen, there are two overall trends. First, using a simpler set of inputs is more robust, especially to neural network architecture perturbations. Second, the overall trend indicates that training on two filters produces results that are better than one-filter training for complex inputs. For simple inputs, the overall trend indicates that two-filter training is relatively similar to one-filter training, as the input distribution shift is reduced due to the use of simpler inputs, although in some cases two-filter training can still produce improvements. In general, applying both methods creates the most robust model. Future work can involve evaluating these approaches on lower-order methods and even more complex flows.

8. Conclusion

Two methods to increase the a posteriori robustness of neural networks are proposed, data augmentation and decreasing the complexity of neural network inputs. Data augmentation involves training with multiple filters in the training data, and two filters have been proposed. The BTF accounts for two-thirds dealiasing as an additional filter, whereas DSCF is non-oscillatory in physical space and uses a single mode approximation to the sharp spectral cutoff filter. A priori, neural networks suffer no performance degradation when trained on one filter as compared with two filters, suggesting that neural network models are able to distinguish between various filters. A posteriori, neural network models trained with two filters are more robust than neural networks trained with one filter, and this trend is confirmed from running neural network models on two high-order solvers. In addition, the neural network trained with three filters is able to be applied to both flow solvers without retraining, and outperforms the neural network trained only with the box filter. Furthermore, neural network models with less complex inputs are more consistent, reducing the distribution shift a posteriori. Training with two different filters often does not reduce the performance of a neural network subgrid stress model either a priori or a posteriori. These trends hold across small neural network architecture perturbations and input normalisations, but are drawn from two high-order solvers and will need to be further evaluated on lower-order solvers. In general, for higher-order solvers, training with data filtered with two different filters and also using less complex inputs to neural networks increases the robustness of neural networks a posteriori.

Acknowledgements

This work used CPU and GPU resources via Bridges2 at the Pittsburgh Supercomputing Center through allocation PHY230025 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by National Science Foundation grants nos. 2138259, 2138286, 2138307, 2137603 and 2138296. The DNS data were provided by the Johns Hopkins Turbulence Database. We thank the anonymous referees whose comments helped improve the paper.

Funding

A.W. is partially supported by NASA Cooperative Agreement number 80NSSC22M0108 and Northrop Grumman, as well as the NDSEG Fellowship.

Declaration of interests

The authors report no conflict of interest.

References

Abekawa, A., Minamoto, Y., Osawa, K., Shimamoto, H. & Tanahashi, M. 2023 Exploration of robust machine learning strategy for subgrid scale stress modeling. Phys. Fluids 35 (1), 015162.CrossRef Google Scholar

Bae, H.J. & Koumoutsakos, P. 2022 Scientific multi-agent reinforcement learning for wall-models of turbulent flows. Nat. Commun. 13 (1), 1443.CrossRef Google Scholar PubMed

Beck, A., Flad, D. & Munz, C.-D. 2019 Deep neural networks for data-driven LES closure models. J. Comput. Phys. 398, 108910.CrossRef Google Scholar

Beck, A. & Kurz, M. 2021 A perspective on machine learning methods in turbulence modelling. Surv. Appl. Maths Mech. 44 (1), e202100002.Google Scholar

Bogacki, P. & Shampine, L.F. 1996 An efficient Runge–Kutta (4,5) pair. Comput. Maths Applics 32 (6), 15–28.CrossRef Google Scholar

Bou-Zeid, E., Meneveau, C. & Parlange, M. 2005 A scale-dependent Lagrangian dynamic model for large eddy simulation of complex turbulent flows. Phys. Fluids 17 (2), 025105.CrossRef Google Scholar

Cheng, Y. et al. 2022 Deep learning for subgrid-scale turbulence modeling in large-eddy simulations of the convective atmospheric boundary layer. J. Adv. Model. Earth Syst. 14 (5), e2021MS002847.CrossRef Google Scholar

Cho, C., Park, J. & Choi, H. 2024 A recursive neural-network-based subgrid-scale model for large eddy simulation: application to homogeneous isotropic turbulence. J. Fluid Mech. 1000, A76.CrossRef Google Scholar

Duraisamy, K. 2021 Perspectives on machine learning-augmented Reynolds-averaged and large eddy simulation models of turbulence. Phys. Rev. Fluids 6 (5), 050504.CrossRef Google Scholar

Ghate, A.S. & Lele, S.K. 2017 Subfilter-scale enrichment of planetary boundary layer large eddy simulation using discrete Fourier–Gabor modes. J. Fluid Mech. 819, 494–539.CrossRef Google Scholar

Guan, Y., Subel, A., Chattopadhyay, A. & Hassanzadeh, P. 2023 Learning physics-constrained subgrid-scale closures in the small-data regime for stable and accurate les. Phys. D: Nonlinear Phenom. 443, 133568.CrossRef Google Scholar

Hu, Z., Subramaniam, A., Kuang, Z., Lin, J., Yu, S., Hannah, W.M., Brenowitz, N.D., Romero, J. & Pritchard, M.S. 2025 Stable machine-learning parameterization of subgrid processes in a comprehensive atmospheric model learned from embedded convection-permitting simulations. arXiv: 2407.00124.CrossRef Google Scholar

Kang, M., Jeon, Y. & You, D. 2023 Neural-network-based mixed subgrid-scale model for turbulent flow. J. Fluid Mech. 962, A38.CrossRef Google Scholar

Kim, H., Shankar, V., Viswanathan, V. & Maulik, R. 2024 Generalizable data-driven turbulence closure modeling on unstructured grids with differentiable physics. Comput. Fluids 318, 107200.CrossRef Google Scholar

Kim, J., Kim, H., Kim, J. & Lee, C. 2022 Deep reinforcement learning for large-eddy simulation modeling in wall-bounded turbulence. Phys. Fluids 34 (10), 105132.CrossRef Google Scholar

Kurz, M., Offenhäuser, P. & Beck, A. 2023 Deep reinforcement learning for turbulence modeling in large eddy simulations. Intl J. Heat Fluid Flow 99, 109094.CrossRef Google Scholar

Li, Y., Perlman, E., Wan, M., Yang, Y., Meneveau, C., Burns, R., Chen, S., Szalay, A. & Eyink, G. 2008 A public turbulence database cluster and applications to study Lagrangian evolution of velocity increments in turbulence. J. Turbul. 9 (31), N31.CrossRef Google Scholar

Ling, J., Kurzawski, A. & Templeton, J. 2016 Reynolds averaged turbulence modelling using deep neural networks with embedded invariance. J. Fluid Mech. 807, 155–166.CrossRef Google Scholar

Maejimam, S. & Kawai, S. 2024 Coarse-grid large-eddy simulation by unsupervised-learning-based sub-grid scale modeling. In AIAA SciTech 2024 Forum, AIAA Paper 2024-1361.Google Scholar

Meneveau, C., Lund, T.S. & Cabot, W.H. 1996 A Lagrangian dynamic subgrid-scale model of turbulence. J. Fluid Mech. 319, 353–385.CrossRef Google Scholar

Meyers, J. & Sagaut, P. 2006 On the model coefficients for the standard and the variational multi-scale Smagorinsky model. J. Fluid Mech. 569, 287–319.CrossRef Google Scholar

Nicoud, F., Toda, H.B., Cabrit, O., Bose, S. & Lee, J. 2011 Using singular values to build a subgrid-scale model for large eddy simulations. Phys. Fluids 23 (8), 085106.CrossRef Google Scholar

Novak, R., Bahri, Y., Abolafia, D.A., Pennington, J. & Sohl-Dickstein, J. 2018 Sensitivity and generalization in neural networks: an empirical study. arXiv: 1802.08760.Google Scholar

Park, J. & Choi, H. 2021 Toward neural-network-based large eddy simulation: application to turbulent channel flow. J. Fluid Mech. 914, A16.CrossRef Google Scholar

Perlman, E., Burns, R., Li, Y. & Meneveau, C. 2007 Data exploration of turbulence simulations using a database cluster. In Supercomputing SC07, ACM. IEEE.CrossRef Google Scholar

Pope, S.B. 2000 Turbulent Flows, 1st edn. Cambridge University.Google Scholar

Sagaut, P. 2006 Large Eddy Simulation for Incompressible Flows: An Introduction, 3rd edn. Springer Science & Business Media.Google Scholar

Sarghini, F., de Felice, G. & Santini, S. 2003 Neural networks based subgrid scale modeling in large eddy simulations. Comput. Fluids 32 (1), 97–108.CrossRef Google Scholar

Sirignano, J. & MacArt, J.F. 2023 Deep learning closure models for large-eddy simulation of flows around bluff bodies. J. Fluid Mech. 966, A26.CrossRef Google Scholar

Sirignano, J., MacArt, J.F. & Freund, J.B. 2020 DPM: a deep learning PDE augmentation method with application to large-eddy simulation. J. Comput. Phys. 423, 109811.CrossRef Google Scholar

Smagorinsky, J. 1963 General circulation experiments with the primitive equations: I. The basic experiment. Mon. Weath. Rev. 91 (3), 99–164.2.3.CO;2>CrossRef Google Scholar

Song, H., Ghate, A.S., Matsuno, K.V., West, J.R., Subramaniam, A. & Lele, S.K. 2024 A robust compact finite difference framework for simulations of compressible turbulent flows. J. Comput. Phys. 519, 113419.CrossRef Google Scholar

Stallcup, E.W. & Dahm, W.J. 2022 Adaptive scale-similar closure for large eddy simulations. Part 1: subgrid stress closure. In AIAA SciTech, AIAA Paper 2022-0595. AIAA.CrossRef Google Scholar

Stoffer, R., van Leeuwen, C.M., Podareanu, D., Codreanu, V., Veerman, M.A., Janssens, M., Hartogensis, O.K. & van Heerwaarden, C.C. 2021 Development of a large-eddy simulation subgrid model based on artificial neural networks: a case study of turbulent channel flow. Geosci. Model Develop. 14 (6), 3769–3788.CrossRef Google Scholar

Wu, A. & Lele, S.K. 2025 a Subgrid stress modelling with multi-dimensional state space sequence models. arXiv: 2511.10910.CrossRef Google Scholar

Wu, A. & Lele, S.K. 2025 b Two neural network unet architecture for subfilter stress modeling. Phys. Rev. Fluids 10 (1), 014601.CrossRef Google Scholar

Xie, C., Wang, J. & Weinan, E. 2020 a Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence. Phys. Rev. Fluids 5 (5), 054606.CrossRef Google Scholar

Xie, C., Yuan, Z. & Wang, J. 2020 b Artificial neural network-based nonlinear algebraic models for large eddy simulation of turbulence. Phys. Fluids 32 (11), 115101.CrossRef Google Scholar

Figure 1. Figure 1 long description.Characterisations of the DSCF in both physical and spectral space. ωc=2πfc$\omega _c = 2 \pi f_c$ corresponds to the filter cutoff, here designated at a 16ΔDNS$16\varDelta _{\textit{DNS}}$ filter width, while ω=2πfn,fn∈[−0.5,0.5]$\omega = 2 \pi f_n, f_n \in [-0.5, 0.5]$ denotes the normalised angular frequency. Two different DSCF widths are shown, one with a 17 and the other with a 33 point support. The DSCF with the smaller filter width has twice the cutoff frequency.

Figure 2. Overall GNN architecture. Note that ‘res. conn.’ denotes residual connections, and these are used in conjunction with skip connections to prevent oversmoothing. There are two layer normalisation layers placed in this architecture. Here ReLU is rectified linear unit.

Figure 3. Graph neural network variant architecture. Note that ‘res. conn.’ denotes residual connections, and these are used in conjunction with skip connections to prevent oversmoothing. There are five layer normalisation layers placed in this architecture, three more than the original. Other than that, the architecture is exactly the same.

Figure 4. Figure 4 long description.Effect of various filters on the ground truth subgrid stress quantities. All colour bars are scaled the same per row: (a) τ11$\tau _{\textit{11}}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations; (b) τ12$\tau _{12}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations; (c) ϵ=−τijS¯ij$\epsilon = -\tau _{\textit{ij}}\bar {S}_{\textit{ij}}$ spatial distribution and PDF. The x and y axis correspond to the coordinates in the physical domain for the spatial distribution visualisations.

Table 1. Neural network configurations; ‘S’ stands for the spectral solver with two-thirds dealiasing (PadeOps), and ‘C’ stands for the sixth-order compact solver (PadeLibs). Normalisation and neural network structure are also given.Table 1 long description.

Table 2. A priori evaluation on the hold-out test set, each value is given as mean (standard deviation). For the RMSE and DRMSE, lower values are better, while higher values are better for the correlation.Table 2 long description.

Table 3. A posteriori evaluation for the Fourier-spectral code for all cases tested.Table 3 long description.

Figure 5. A posteriori one- and two-filter comparisons with traditional models for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Figure 6. Figure 6 long description.A posteriori multiple filter comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation, but trained with a variety of filter combinations ranging from one to three different filter shapes.

Figure 7. A posteriori ALN comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the variant ALN architecture with the global normalisation.

Figure 8. Figure 8 long description.A posteriori simple input comparisons for the Fourier-spectral code with two-thirds dealiasing. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. Both the original and ALN GNN architectures are tested with ‘simple’ inputs using both local and global normalisation.

Table 4. A posteriori evaluation for the sixth-order compact code for all cases tested.Table 4 long description.

Figure 9. Figure 9 long description.A posteriori one- and two-filter comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Figure 10. A posteriori multiple-filter comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation, but trained with a variety of filter combinations ranging from one to three different filter shapes.

Figure 11. Figure 11 long description.A posteriori ALN comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the ALN architecture with the global normalisation.

Figure 12. A posteriori simple input comparisons for the sixth-order compact code. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. Both the original and ALN GNN architectures are used with the global and local normalisation for the ‘simple’ inputs.

Figure 13. Figure 13 long description.A posteriori one- versus two-filter comparison with the Sigma model as a reference for a coarser grid spacing corresponding to ΔLES=32ΔDNS$\varDelta _{\textit{LES}} = 32 \varDelta _{\textit{DNS}}$. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation.

Figure 14. A posteriori simple input comparisons with the Sigma model as a reference on the coarser grids. Spectra, longitudinal flatness and the PDF of ωi¯Sij¯ωj¯$\overline {\omega _i}\overline {S_{\textit{ij}}}\overline {\omega _{\!j}}$ are plotted. The GNN architecture used corresponds to the original architecture with the global normalisation applied to the ‘simple’ inputs.

Figure 15. Figure 15 long description.A posteriori streamwise velocity spectra and mean profiles for channel flow, comparing one- versus two-filter GNN models with traditional models. The GNN architecture used corresponds to the original architecture with the global normalisation.

Figure 16. A posteriori one- versus two-filter Reynolds stress comparisons with traditional models for channel flow. Here, z/H=0$z/H = 0$ corresponds to the wall while z/H=1$z/H=1$ corresponds to the channel centreline. The GNN architecture used corresponds to the original architecture with the global normalisation.

Figure 17. Figure 17 long description.A posteriori spectra and mean profile comparison of different neural network perturbations. Both the original and ALN GNN architectures are used, with global normalisation applied to the ‘complex’ inputs.

Figure 18. A posteriori Reynolds stress comparison of different neural network perturbations. Here, z/H=0$z/H = 0$ corresponds to the wall while z/H=1$z/H=1$ corresponds to the channel centreline. Both the original and ALN GNN architectures are used, with global normalisation applied to the ‘complex’ inputs.

Figure 19. A posteriori simple inputs streamwise velocity spectra and mean profiles. The GNN original architecture is used, with global normalisation applied to the ‘simple’ inputs.

Figure 20. A posteriori simple inputs Reynolds stress profiles for channel flow. Here, z/H=0$z/H = 0$ corresponds to the wall while z/H=1$z/H=1$ corresponds to the channel centreline. The GNN original architecture is used, with global normalisation applied to the ‘simple’ inputs.

Article contents

Addressing a posteriori performance degradation in neural network subgrid stress models

Abstract

JFM classification

Information

1. Introduction

2. Numerical details

2.1. Large eddy simulation

2.2. A posteriori flow solvers

3. Filtering and neural network details

3.1. Discrete spatial cosine filter: a localised low-pass filter

3.1.1 Two-thirds dealiasing as an additional filter

3.1.2. Neural network details

3.2. Training data

3.2.1. Training configuration

3.2.2. Effect of filtering on ground truth

3.3. All neural network configurations

4. A priori analysis

4.1. One-filter versus two-filter results

5. A posteriori analysis

5.1. Fourier-spectral with two-thirds dealiasing (PadeOps)

5.1.1. Effect of multiple filters

5.1.2. Neural network perturbation

5.1.3. Using simpler inputs

5.2. Sixth-order compact (PadeLibs)

5.2.1. Effect of multiple filters

5.2.2. Neural network perturbation

5.2.3. Using simpler inputs

6. Additional a posteriori evaluations

6.1. Evaluation on a different grid width

6.2. Evaluation on channel flow

6.2.1. Two-filter complex input analysis

6.2.2. Two-filter complex input analysis ALN

6.2.3. Simple inputs

7. Combined a posteriori analysis

8. Conclusion

Acknowledgements

Funding

Declaration of interests

References

Save article to Kindle

Save article to Dropbox

Save article to Google Drive

Reply to: Submit a response

Your details

You have entered the maximum number of contributors

Conflicting interests