2. Learning turbulence models

In this paper, we study neural networks for turbulence modelling in incompressible fluids. These flows are governed by the Navier–Stokes equations

(2.1) \begin{equation} \left.\begin{array}{c@{}} \dfrac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u} ={-}\boldsymbol{\nabla} p + \dfrac{1}{{Re}}\nabla^{2}\boldsymbol{u} + \boldsymbol{f} ,\\ \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{array}\right\} \end{equation}

where $\boldsymbol {u}= [u \ v]^\textrm {T}$, $p$ and $Re$ are the velocity field, pressure field and Reynolds number respectively. The term $\boldsymbol {f} = [f_x \ f_y]^\textrm {T}$ represents an external force on the fluid. In the context of turbulent flows, an accurate solution to these equations entails either resolving and numerically simulating all turbulent scales, or modelling the turbulence physics through an approximative model.

Our aim is to develop a method that enhances fluid simulations by means of a machine learning model. In particular, we aim to improve the handling of fine temporal and spatial turbulence scales that are potentially under-resolved, such that the influence of these scales on the larger resolved motions needs to be modelled. The function that approximates these effects is solely based on low-resolution data and is herein parameterised by a CNN. The network is then trained to correct a low-resolution numerical solution during the simulation, such that the results coincide with a downsampled high-resolution dataset. Within this hybrid approach, the turbulence model directly interacts with the numerical solver at training and at inference time. To achieve this objective, we utilise differentiable solvers, i.e. solvers which provide derivatives with respect to their output state. Such solvers can be seen as part of the differentiable programming methodology in deep learning, which is equivalent to employing the adjoint method from classical optimisation (Giles et al. Reference Giles, Duta, Muller and Pierce2003) in the context of neural networks. The differentiability of the solver enables the propagation of optimisation gradients through multiple solver steps and neural network evaluations.

2.1. Differentiable PISO solver

Our differentiable solver is based on the semi-implicit pressure-implicit with splitting of operators (PISO) scheme introduced by Issa (Reference Issa1986), which has been used for a wide range of flow scenarios (Kim & Benson Reference Kim and Benson1992; Barton Reference Barton1998). Each second-order time integration step is split into an implicit predictor step solving the discretised momentum equation, followed by two corrector steps that ensure the incompressibility of the numerical solution to the velocity field. The Navier–Stokes equations are discretised using the finite-volume method, while all cell fluxes are computed to second-order accuracy.

The solver is implemented on the basis of TensorFlow (TF) (Abadi Reference Abadi2016), which facilitates parallel execution of linear algebra operations on the graphics processing unit (GPU), as well as the differentiability of said operations. Additional functions exceeding the scope of TF are written as custom operations and implemented using compute unified device architecture (CUDA) programming. This approach allows us to seamlessly integrate initially unsupported features such as sparse matrix operations in the TF graph. More details about the solver can be found in Appendix A, where the solver equations are listed in Appendix A.1, implementation details in Appendix A.2 and a verification is conducted in Appendix A.3. Figure 1 gives a brief overview of the solver procedure.

Figure 1. Solver procedure of the PISO scheme and its interaction with the convolutional neural network; data at time $t_n$ are taken from the DNS dataset and processed by the downsampling operation q that yields downsampled representations of the input fields, before entering the differentiable solver; the solver unrollment performs $m$ steps, each of which is corrected by the CNN, and is equivalent to $\tau$ high-resolution steps; the optimisation loss takes all resulting (intermediate) timesteps.

In the following, we will denote a PISO solver step $\mathcal {S}$ as

(2.2)\begin{equation} (\boldsymbol{u}_{n+1}, p_{n+1}) = \mathcal{S}(\boldsymbol{u}_n, p_n, \boldsymbol{f}_n), \end{equation}

where $\boldsymbol {u}_{n}$, $p_{n}$ and $\boldsymbol {f}_n$ represent discretised velocity, pressure and forcing fields at time $t_n$.

2.3. Unrolling timesteps for training

Our method combines the numerical solver introduced in § 2.1 with the modelling capabilities of CNNs as outlined in § 2.2. As also illustrated in figure 1, the resulting data-driven training algorithm works based on a dataset $(\boldsymbol {u}(t_n), p(t_n))$ consisting of high-resolution ($N_x\times N_y$) velocity fields $\boldsymbol {u}(t_n)\in \mathbb {R}^{N_x\times N_y \times 2}$ and corresponding pressure fields $p(t_n)\in \mathbb {R}^{N_x\times N_y}$ for the discrete time $t_n$. In order to use these DNS data for training under-resolved simulations on different grid resolutions, we define a downsampling procedure $q(\boldsymbol {u}, p):\mathbb {R}^{ N_x\times N_y \times 3}\xrightarrow {}\mathbb {R}^{\tilde N_x \times \tilde N_y\times 3}$ , that takes samples from the dataset and outputs the data $(\tilde {\boldsymbol {u}}_n,\tilde {p}_n)$ at a lower target resolution ($\tilde N_x\times \tilde N_y$) via bilinear interpolation. This interpolation provides a simple method of acquiring data at the shifted cell locations of different discretisations. It can be seen as a repeated linear interpolation to take care of two spatial dimensions. The resampling of DNS data is used to generate input and target frames of an optimisation step. For the sake of simplicity, we will denote a downsampled member of the dataset consisting of velocities and pressure as $\tilde {q}_n=q(\boldsymbol {u}(t_n), p(t_n))$. Similarly, we will write $\tilde {\boldsymbol {f}}_n= \boldsymbol {f}_{CNN}(\tilde {\boldsymbol {u}}_{n},\boldsymbol {\nabla } \tilde {p}_{n}|\theta )$. Note that the network operates solely on low-resolution data and introduces a corrective forcing to the low-resolution simulation, with the goal of reproducing the behaviour of a DNS. We formulate the training objective as

(2.3)\begin{equation} \min_\theta(\mathcal{L}(\tilde{q}_{n+\tau}, \mathcal{S}_\tau(\tilde{q}_n,\tilde{\boldsymbol{f}}_n))), \end{equation}

for a loss function $\mathcal {L}$ that satisfies $\mathcal {L}(x,y)\xrightarrow {}0$ for $x\approx y$. By this formulation, the network takes a downsampled DNS snapshot and should output a forcing which makes the flow fields after a low-resolution solver step closely resemble the next downsampled frame. The temporal increment $\tau$ between these subsequent frames is set to match the timesteps in the low-resolution solver $\mathcal {S}$, which in turn are tuned to maintain Courant numbers identical to the DNS.

Um et al. (Reference Um, Brand, Fei, Holl and Thuerey2020) showed that similar training tasks benefit from unrolling multiple temporal integration steps in the optimisation loop. The optimisation can then account for longer-term effects of the network output on the temporal evolution of the solution, increasing accuracy and stability in the process. We utilise the same technique and find it to be critical for the long-term stability of turbulence models. Our notation from equations (2.2) and (2.3) is extended to generalise the formulation towards multiple subsequent snapshots. When training a model through $m$ unrolled steps, the optimisation objective becomes

(2.4)\begin{equation} \min_\theta\left(\sum_{s=0}^{m}\mathcal{L}(\tilde{q}_{n+s\tau} ,\mathcal{S}_\tau^{s}(\tilde{q}_n,\tilde{\boldsymbol{f}}_n))\right), \end{equation}

where $\mathcal {S}^{s}$ denotes the successive execution of $s$ solver steps including network updates, starting with the initial fields $q_i$. By this formulation the optimisation works towards matching not only the final, but also all intermediate frames. Refer to Appendix A.1 for a detailed explanation of this approach, including equations for optimisation and loss differentiation.

2.4. Loss functions

As introduced in (2.3), the training of deep CNNs is an optimisation of its parameters. The loss function $\mathcal {L}$ serves as the optimisation objective and thus has to assess the quality of the network output. Since our approach targets the reproduction of DNS-like behaviour on a coarse gird, the chosen loss function should consequently aim to minimise the distance between the state of a modelled coarse simulation and the DNS. In this context, a natural choice is the $\mathcal {L}_2$ loss on the $s$th unrolled solver step

(2.5)\begin{equation} \mathcal{L}_2 = \sqrt{(\tilde{\boldsymbol{u}}_{s}- q(\boldsymbol{u}_{s\tau}))\boldsymbol{\cdot}(\tilde{\boldsymbol{u}}_{s}- q(\boldsymbol{u}_{s\tau}))}, \end{equation}

since this formulation drives the optimisation towards resembling a desired outcome. Therefore, the $\mathcal {L}_2$ loss trains the network to directly reproduce the downsampled high-resolution fields, and the perfect reproduction from an ideal model gives $\mathcal {L}_2=0$. Since the differentiable solver allows us to unroll multiple simulation frames, we apply this loss formulation across a medium-term time horizon and thus also optimise towards multi-step effects. By repeatedly taking frames from a large DNS dataset in a stochastic sampling process, a range of downsampled instances are fed to the training procedure. While the DNS dataset captures all turbulence statistics, they are typically lost in an individual training iteration. This is due to the fact that training mini-batches do not generally include sufficient samples to represent converged statistics, and no specific method is used to satisfy this criterion. This means that data in one training iteration solely carry instantaneous information. Only the repeated stochastic sampling from the dataset lets the network recover awareness of the underlying turbulence statistics. The repeated matching of instantaneous DNS behaviour thus encodes the turbulence statistics in the training procedure. While the $\mathcal {L}_2$ loss described in (2.5) has its global minimum when the DNS behaviour is perfectly reproduced, in practice, we find that it can neglect the time evolution of certain fine scale, low amplitude features of the solutions. This property of the $\mathcal {L}_2$ loss is not unique to turbulence modelling and has previously been observed in machine learning in other scientific fields such as computer vision (Yu et al. Reference Yu, Fernando, Hartley and Porikli2018). To alleviate these shortcomings, we include additional loss formulations, which alter the loss landscape to avoid these local minima.

We define a spectral energy loss $\mathcal {L}_{E}$, designed to improve the accuracy of the learned model on fine spatial scales. It is formulated as the log-spectral distance of the spectral kinetic energies at the $s$th step

(2.6)\begin{equation} \mathcal{L}_E = \sqrt{\int_k \log\left(\frac{\tilde{E}_{s}(k)}{E_{s \tau}^{q}(k)}\right)^{2} \, {{\rm d}}k}, \end{equation}

where $\tilde {E}_s(k)$ denotes the spectral kinetic energy of the low-resolution velocity field at wavenumber $k$, and $E^{q}_{s\tau }$ represents the same quantity for the downsampled DNS data. In practice, this loss formulation seeks to equalise the kinetic energy in the velocity field for each discrete wavenumber. The $\log$ rescaling of the two respective spectra regularises the relative influence of different spatial scales. This energy loss elevates the relative importance of fine scale features.

Our final aim is to train a model that can be applied to a standalone forward simulation. The result of a neural network-modelled low-resolution simulation step should therefore transfer all essential turbulence information, such that the same model can once again be applied in the subsequent step. The premises of modelling the unresolved behaviour are found in the conservation equation for the implicitly filtered low-resolution kinetic energy in tensor notation

(2.7)\begin{equation} \frac{\partial \widetilde{E_f}}{\partial t} + \tilde{u}_i \frac{\partial \widetilde{E_f}}{\partial x_i}+\frac{\partial}{\partial x_j} \tilde{u}_i \left[ \delta_{ij}\tilde{p}+ \tau_{ij}^{r} - \frac{2}{{Re}}\tilde{{\mathsf{S}}}_{ij} \right] ={-}\epsilon_f-\mathcal{P}^{r}, \end{equation}

where $\tilde {E}_f$ denotes the kinetic energy of the filtered velocity field, $\tau _{ij}^{r}$ represents the SGS stress tensor, ${\tilde {{\mathsf{S}}}}_{ij}= \frac {1}{2}({\partial \tilde {u}_i}/{\partial x_j} + {\partial \tilde {u}_j}/{\partial x_i})$ is the resolved rate of strain and $\epsilon _f$ and $\mathcal {P}^{r}$ are sink and source terms for the filtered kinetic energy. These terms are defined as $\epsilon _f = ({2}/{Re}){\tilde {{\mathsf{S}}}}_{ij}{\tilde {{\mathsf{S}}}}_{ij}$ and $\mathcal {P}^{r}=-\tau _{ij}^{r}{\tilde {{\mathsf{S}}}}_{ij}$. The viscous sink $\epsilon _f$ represents the dissipation of kinetic energy due to molecular viscous stresses at grid-resolved scales. In hybrid methods, this viscous dissipation at grid level $\epsilon _f$ is fully captured by the numerical solver. On the contrary, the source term $\mathcal {P}^{r}$ representing the energy transfer from resolved scales to residual motions cannot be computed, because the SGS stresses $\tau _{ij}^{r}$ are unknown. One part of the modelling objective is to estimate these unresolved stresses and the interaction of filtered and SGS motions. Since the energy transfer between these scales $\mathcal {P}^{r}$ depends on the filtered rate of strain ${\tilde{{\mathsf{S}}}}_{ij}$, special emphasis is required to accurately reproduce the filtered rate of strain tensor. This motivates the following rate of strain loss at the $s$th unrolled solver step

(2.8)\begin{equation} \mathcal{L}_{{S}}= \sum_{i,j}|\tilde{{{\mathsf{S}}}}_{ij, s} - {{\mathsf{S}}}^{q}_{ij, s\tau}| , \end{equation}

where ${{\mathsf{S}}}^{q}_{ij, s}$ denotes the rate of strain of the downsampled high-resolution velocity field. This loss term ensures that the output of a hybrid solver step carries the information necessary to infer an accurate forcing in the subsequent step.

While our models primarily focus on influences of small-scale motions on the large-scale resolved quantities, we now draw attention to the largest scale, the mean flow. To account for the mean flow at training time, an additional loss contribution is constructed to match the multi-step statistics and written as

(2.9)\begin{equation} \mathcal{L}_{MS}= \| \langle \boldsymbol{u}_s \rangle_{s=0}^{m} - \langle \tilde{\boldsymbol{u}}_{s\tau}\rangle_{s=0}^{m}\|_1, \end{equation}

where $\langle \rangle _{s=0}^{m}$ denotes an averaging over the $m$ unrolled training steps with iterator $s$. This notation resembles Reynolds averaging, albeit being focused on the shorter time horizon unrolled during training. Matching the averaged quantities is essential to achieving long-term accuracy of the modelled simulations for statistically steady simulations, but lacks physical meaning in transient cases. Therefore, this loss contribution is solely applied to statistically steady simulations. In this case, the rolling average $\langle \rangle _{s=0}^{m}$ approaches the steady mean flow for increasing values of $m$.

Our combined turbulence loss formulation as used in the network optimisations additively combines the aforementioned terms as

(2.10)\begin{equation} \mathcal{L}_{T} = \lambda_2\mathcal{L}_2 + \lambda_E\mathcal{L}_E + \lambda_{S}\mathcal{L}_{\boldsymbol{S}} + \lambda_{MS}\mathcal{L}_{MS}, \end{equation}

where $\lambda$ denotes the respective loss factor. Their exact values are mentioned in the flow scenario specific chapters. Note that these loss terms, similar to the temporal unrolling, do not influence the architecture or computational performance of the trained neural network at inference time. They only exist at training time to guide the network to an improved state with respect to its trainable parameters. In the following main sections of this paper, we use three different turbulence scenarios to study effects of the number of unrolled steps and the whole turbulence loss $\mathcal {L}_{T}$. An ablation on the individual components of $\mathcal {L}_{T}$ is provided in Appendix F. We start with employing the training strategy on isotropic decaying turbulence.

3. Two-dimensional isotropic decaying turbulence

Isotropic decaying turbulence in two dimensions provides an idealised flow scenario (Lilly Reference Lilly1971), and is frequently used for evaluating model performance (San Reference San2014; Maulik et al. Reference Maulik, San, Rasheed and Vedula2019; Kochkov et al. Reference Kochkov, Smith, Alieva, Wang, Brenner and Hoyer2021). It is characterised by a large number of vortices that merge at the large spatial scales whilst the small scales decay in intensity over time. We use this flow configuration to explore and evaluate the relevant parameters, most notably the number of unrolled simulation steps as well as the effects of loss formulations.

In order to generate training data, we ran a simulation on a square domain with periodic boundaries in both spatial directions. The initial velocity and pressure fields were generated using the initialisation procedure by San & Staples (Reference San and Staples2012). The Reynolds numbers are calculated as $Re = (\hat {e} \hat {l})/\nu$, with the kinetic energy $\hat {e}=(\langle u_iu_i\rangle )^{1/2}$ and the integral length scale $\hat {l} = \hat {u}/\hat {\omega }$ and $\hat {\omega }=(\langle \omega _i\omega _i\rangle )^{1/2}$. The Reynolds number of this initialisation was $Re=126$. The simulation was run for a duration of $T=10^{4}\Delta t_{{DNS}}= 100\hat {t}$, where the integral time scale is calculated as $\hat {t}=1/\hat {\omega }$ at the initial state. During the simulation, the backscatter effect transfers turbulence energy to the larger scales, which increases the Reynolds number (Kraichnan Reference Kraichnan1967; Chasnov Reference Chasnov1997). In our dataset, the final Reynolds number was $Re=296$. Note that despite this change in Reynolds number, the turbulence kinetic energy is still decreasing and the flow velocities will decay to zero.

Our aim is to estimate the effects of small-scale turbulent features on a coarser grid based on fully resolved simulation data. Consequently, the dataset should consist of a fully resolved DNS and satisfy the resolution requirements. In this case the square domain $(L_x,L_y)=(2{\rm \pi},2{\rm \pi} )$ was discretised by $(N_x,N_y) = (1024,1024)$ grid cells and the simulation evolved with a timestep satisfying CFL $=0.3$.

We trained a series of models on downsampled data, where spatial and temporal resolutions were decreased by a factor of $8$, resulting in an effective simulation size reduction of $8^{3}=512$. Our best performing model was trained through $30$ unrolled simulation steps. This is equivalent to $1.96\hat t$ for the initial simulation state. Due to the decaying nature of this test case, the integral time scale increases over the course of the simulation, while the number of unrolled timesteps is kept constant. As a consequence, the longest unrollments of 30 steps cover a temporal horizon similar to the integral time scale. All our simulation set-ups will study unrollment horizons ranging up to the relevant integral time scale, and best results are achieved when the unrollment approaches the integral time scale. For training the present set-up, the loss factors from equation (2.10) were chosen as $(\lambda _2 , \lambda _E, \lambda _{S}, \lambda _{MS})= (10, 5\times 10^{-2}, 1\times 10^{-5}, 0)$.

To evaluate the influence of the choice of loss function and the number of unrolled steps, several alternative models were evaluated. Additionally, we trained a model with a traditional supervised approach. In this setting, the differentiable solver is not used, and the training is performed purely on the basis of the training dataset. In this case, the corrective forcing is added after a solver step is computed. The optimisation becomes

(3.1)\begin{equation} \min_\theta (\mathcal{L}(q_{n+\tau},\boldsymbol{f}_{CNN} (\mathcal{S}_\tau(q_n))). \end{equation}

The equations for the supervised training approach are detailed in Appendix A.1. Furthermore, a LES with the standard Smagorinsky model was included in the comparison. A parameter study targeting the Smagorinsky coefficient revealed that a value of $C_s=0.008$ handles the physical behaviour of our set-up best. See Appendix D for details. An overview of all models and their parameters is given in table 1.

Table 1. Training details for models trained on the isotropic turbulence case, wtih NN representing various versions of our neural network models; mean squared error (MSE) evaluated at $t_1 = 64\Delta t\approx 5\hat {t}$ and $t_2 = 512\Delta t\approx 40\hat {t}$.

After training, a forward simulation was run for comparison with a test dataset. For the test data, an entirely different, randomly generated initialisation was used, resulting in a velocity field different from the simulation used for training. The test simulations were advanced for $1000\Delta t = 80\hat t$.

Note that the temporal advancement of the forward simulations greatly surpasses the unrolled training horizon, which leads to instabilities with the supervised and 1-step model, and ultimately to the divergence of their simulations. Consequently, we conclude that more unrolled steps are critical for the applicability of the learned models and do not include the 1-step model in further evaluations. While an unrollment of multiple steps also improves the accuracy of supervised models, these models nevertheless fall short of their differentiable counterparts, as seen in a deeper study in Appendix E.

We provide visual comparisons of vorticity snapshots in figure 2, where our method's improvements become apparent. The network-modelled simulations produce a highly accurate evolution of vorticity centres, and comparable performance cannot be achieved without a model for the same resolution. We also investigate the resolved turbulence kinetic energy spectra in figure 3. Whilst the no-model simulation overshoots the DNS energy at its smallest resolved scales, the learned model simulations perform better and match the desired target spectrum. Figure 4 shows temporal evolutions of the domain-wide-resolved turbulence energy and the domain-wide-resolved turbulence dissipation rate. The turbulence energy is evaluated according to $E(t)= \int u_i'(t) u_i'(t) \, {\textrm {d}} \varOmega$, where $u_i'$ is the turbulent fluctuation. We calculate the turbulence dissipation as $\epsilon (t)=\int \langle \mu ({\partial u'_i\partial u_i'}/{\partial x_j \partial x_j}) \rangle \, {\textrm {d}} \varOmega$. Simulations with our CNN models strongly agree with the downsampled DNS.

Figure 2. Vorticity visualisations of DNS, no-model, LES, and learned model simulations at $t = (350,700) \Delta t$ on the test dataset, zoomed-in version below.

Figure 3. Resolved turbulence kinetic energy spectra of the downsampled DNS, no-model, LES and learned model simulations; the learned 30-step model matches the energy distribution of downsampled DNS data; the vertical line represents the Nyquist wavenumber of the low-resolution grid.

Figure 4. Comparison of DNS, no-model, LES and learned model simulations with respect to resolved turbulence kinetic energy over time (a); and turbulence dissipation rate over time (b).

All remaining learned models stay in close proximity to the desired high-resolution evolutions, whereas the LES-modelled and no-model simulations significantly deviate from the target. Overall, the neural network models trained with more unrolled steps outperformed others, while the turbulence loss formulation $\mathcal {L}_{T}$ also had a positive effect.

In particular, the backscatter effect is crucial for simulations of decaying turbulence (Kraichnan Reference Kraichnan1967; Smith, Chasnov & Waleffe Reference Smith, Chasnov and Waleffe1996). The CNN adequately dampens the finest scales, as seen in the high wavenumber section of the energy spectrum (figure 3), it also successfully boosts larger scale motions. In contrast, the no-model simulation lacks dampening at the finest scales and cannot reproduce the backscatter effect on the larger ones. On the other hand, the dissipative nature of the Smagorinsky model used in the LES leads to under-sized spectral energies across all scales. Especially the spectral energies of no-model and LES around wavenumber $k=10$ show large deviations form the ground truth, while the CNN model accurately reproduces its behaviour. These large turbulent scales are the most relevant to the resolved turbulence energy and dissipation statistics, which is reflected in figure 4. Herein, the neural network models maintain the best approximations, and high numbers of unrolled steps show the best performance at long simulation horizons. The higher total energy of the neural network modelled simulations can be attributed to the work done by the network forcing, which is visualised together with the SGS stress tensor work from the LES simulation as well as its SGS energy in figure 5. This analysis reveals that the neural networks do more work on the system as the LES model does, which explains the higher and more accurate turbulence energy in figure 4 and the spectral energy behaviour at large scales in figure 3.

Figure 5. NN-model work on the flow field, work by the LES model and the estimated SGS energies from LES.

4. Temporally developing planar mixing layers

Next, we apply our method to the simulation of two-dimensional planar mixing layers. Due to their relevance to applications such as chemical mixing and combustion, mixing layers have been the focus of theoretical and numerical studies in the fluid-mechanics community. These studies have brought forth a large set and good understanding of a posteriori evaluations, like the Reynolds-averaged turbulent statistics or the vorticity and momentum thickness. Herein, we use these evaluations to assess the accuracy of our learned models with respect to metrics that are not directly part of the learning targets.

Temporally evolving planar mixing layers are the simplest numerical representation of a process driven by the Kelvin–Helmholtz instability in the shear layer. They are sufficiently defined by the Reynolds number, domain sizes, boundary conditions and an initial condition. Aside from the shear layer represented by a $\tanh$-profile, the initial flow fields feature an oscillatory disturbance that triggers the instability, leading to the roll up of the shear layer. This has been investigated by theoretical studies involving linear stability analysis (Michalke Reference Michalke1964) and numerical simulation (Rogers & Moser Reference Rogers and Moser1994). Our set-up is based on the work by Michalke (Reference Michalke1964), who studied the stability of the shear layer and proposed initialisations that lead to shear layer roll up. As initial condition, we add randomised modes to the mean profile, resulting in the streamfunction

(4.1)\begin{equation} \varPsi(x,y) = y + \tfrac{1}{2}\ln(1+{\rm e}^{{-}4y})+a((\alpha y)^{2}+1)\, {\rm e}^{-(\alpha y)^{2}}\cos(\omega_\varPsi x), \end{equation}

where $a$ is the amplitude of the perturbation, $\alpha$ parameterises the decay of the perturbation in $y$-direction and $\omega _\varPsi$ represents the perturbation frequency. The initial flow field can then be calculated by

(4.2a,b)\begin{equation} u(x,y) = \frac{\partial \varPsi}{\partial y}, \quad v(x,y) ={-} \frac{\partial \varPsi}{\partial x}. \end{equation}

At the initial state this results in a velocity step $\Delta U= {U_2}-{U_1}=1$ and a vorticity thickness of $\delta _\omega ={\Delta U}/({\partial U }/{\partial y}|_{{max}})= 1$, where velocities marked as $U$ represent mean-stream quantities. Thus, $U_2$ and $U_1$ are the fast and slow mean velocities of the shear layer. The computational domain of size $(L_x,L_y)=(40{\rm \pi},20{\rm \pi} )$ is discretised by $(N_x,N_y)=(1024,512)$ grid cells for the high-resolution dataset generation. The streamwise boundaries are periodic, while the spanwise boundaries in the $y$-direction are set to a free-slip boundary, where ${\partial u}/{\partial y} |_{\varOmega _y}=0$, $v|_{\varOmega _y}=0$ and $p|_{\varOmega _y}=0$. The Reynolds number based on the unperturbed mean profile and the vorticity thickness is calculated to be $Re={\Delta U \delta _\omega }/{\nu }=250$ for all randomised initialisations. The simulations are run for $T=420=12\ 000 \Delta t_{DNS}$. Our dataset consists of three simulations based on different initialisations. Their perturbation details are found in table 2. Two of these simulations were used as training datasets, while all of our evaluation is performed on the remaining one as the extrapolation test dataset.

Table 2. Perturbation details for initial conditions of temporal mixing layer training and test datasets.

Following the approach in § 3, the model training uses a $8\times$ downscaling in space and time. The loss composition was set to $(\lambda _2 , \ \lambda _E, \ \lambda _{S}, \ \lambda _{MS})= (100, \ 2, \ 5\times 10^{-2}, \ 0)$. We used the same CNN architecture as introduced earlier, although due to the difference in boundary conditions a different padding procedure was chosen (see Appendix B). To illustrate the impact of the turbulence loss $\mathcal {L}_{T}$ and an unrolling of $60$ numerical steps, we compare with several variants with reduced loss formulations and fewer unrolling steps. The maximum number of $60$ unrolled steps corresponds to $16 t_{\delta _\theta }$ integral time scales computed on the momentum thickness as $t_{\delta _\theta } = \delta _\theta /\Delta U$. With the shear layer growing, the momentum thickness increases $7$-fold, which decreases the number of integral time scales to $2$ for 60 steps of unrollment. Table 3 shows details of the model parameterisations. To avoid instabilities in the gradient calculation that could ultimately lead to unstable training, we split the back-propagation into subranges for the $60$-step model. This method stabilises an otherwise unstable training of the $60$-step model, and a split into $30$-step long back-propagation subranges performs best. Such a model is added to the present evaluations as $\textrm {NN}_{60,\mathcal {L}_{T}}$. Detailed results regarding the back-propagation subranges are discussed in § 6.

Table 3. Model details for unrollment study; MSE with respect to DNS from test data at $t_e = 512\Delta t$.

The trained models were compared with a downsampled DNS and a no-model simulation, all sharing the same starting frame from the test dataset. This test dataset changes the initial condition, where different perturbation frequencies and amplitudes result in a variation in vortex roll-up and vortex merging behaviour of the mixing layer. The resulting numerical solutions were compared at three different evolution times $t=[256 \ 640 \ 1024]\Delta t$. Figure 6 shows the vorticity heat map of the solutions. Qualitatively, the simulations corrected by the CNN exhibit close visual proximity to the DNS by boosting peaks in vorticity where applicable, and additionally achieve a dampening of spurious oscillations.

Figure 6. Vorticity visualisations of DNS, no-model and learned model simulations at $t=(256,640,1024)\Delta t$ on the test dataset.

These observations are matched by corresponding statistical evaluations. The statistics are obtained by averaging the simulation snapshots along their streamwise axis and the resulting turbulence fluctuations were processed for each evaluation time. Figure 7 shows that all $\mathcal {L}_{T}$-models closely approximate the DNS reference with respect to their distribution of resolved turbulence kinetic energy and Reynolds stresses along the cross-section, while the no-model simulation clearly deviates. Note that the mixing process causes a transfer of momentum from fast to slow moving sections through the effects of turbulent fluctuations. The shear layer growth is thus dominated by turbulent diffusion. Consequently, accurate estimates of the turbulent fluctuations are necessary for the correct evolution of the mixing layer. These fluctuations are most visible in the Reynolds stresses $u'v'$, and an accurate estimation is an indicator for well-modelled turbulent momentum diffusion. The evaluations also reveal that unrolling more timesteps during training gains additional performance improvements. These effects are most visible when comparing the 10-step and 60-step model in a long temporal evolution, as seen in the Reynolds stresses in figure 7. The evaluation of resolved turbulence kinetic energies shows that the models correct for the numerical dissipation of turbulent fluctuations, while, in contrast, there is an underestimation of kinetic energy in the no-model simulation. While longer unrollments generally yield better accuracy, it is also clear that $30$ steps come close to saturating the model performance in this particular flow scenario. With the integral time scales mentioned earlier, it becomes clear that 30 simulation steps capture one integral time scale of the final simulation phase, i.e. the phase of the decaying simulation that exhibits the longest time scales. One can conclude that an unrollment of one time scale is largely sufficient, and further improvements of unrolling 2 time scales with $60$ steps are only minor.

Figure 7. Comparison of DNS, no-model and learned model simulations with respect to resolved turbulence kinetic energy (a) and Reynolds stresses (b).

The resolved turbulence kinetic energy spectra are evaluated to assess the spatial scales at which the corrective models are most active. The spectral analysis at the centreline is visualised in figure 8, whilst the kinetic energy obtained from fluctuations across the cross-section with respect to streamwise averages is shown in figure 9. These plots allow two main observations: firstly, the deviation of kinetic energy mostly originates from medium-sized spatial scales, which are dissipated by the no-model simulation but are accurately reconstructed by the neural network trained with $\mathcal {L}_T$. This effect is connected to the dampening of vorticity peaks in the snapshots in figure 6. Secondly, the fine-scale spectral energy of the no-model simulation has an amplitude similar to the DNS over long temporal horizons (figure 9). This can be attributed to numerical oscillations rather than physical behaviour. These numerical oscillations, as also seen in the snapshots in figure 6, exist for the no-model simulation but are missing in the $\mathcal {L}_T$-modelled simulations. Training a model without the additional loss terms in $\mathcal {L}_T$ from (2.10), i.e. only with the $\mathcal {L}_2$ from (2.5), yields a model that is inaccurate and results in unphysical oscillations. It does not reproduce the vorticity centres, and is also unstable over long temporal horizons. Herein, non-physical oscillations are introduced, which also show up in the cross-sectional spectral energies and vorticity visualisations. We thus conclude that best performance can be achieved with a network trained with $\mathcal {L}_T$, which learns to dampen numerical oscillations and reproduces physical fluctuations across all spatial scales.

Figure 8. Centreline kinetic energy spectra for the downsampled DNS, no-model and learned model simulations.

Figure 9. Cross-sectional kinetic energy spectra of the downsampled DNS, no-model and learned model simulations.

It is worth noting that our method is capable of enhancing an under-resolved simulation across a wide range of turbulent motions. The vortex size in the validation simulation ranges from $7\delta _{\omega _0}$ at the starting frame to $60\delta _{\omega _0}$ after evolving for $1200 \Delta t$. This timespan encompasses two vortex merging events, both of which cannot be accurately reproduced with a no-model, or a $\mathcal {L}_2$-model simulation, but are captured by the $\mathcal {L}_T$-trained network models. This is shown in the comparison of the momentum thicknesses over time in figure 10. The reproduction of turbulence statistics (figure 7) yields, in the long term, an accurate turbulent diffusion of momentum and mixing layer growth for the models trained with $\mathcal {L}_T$. On the contrary, the $\mathcal {L}_2$ model fails to reproduce the vortex cores and deviates with respect to the momentum thickness for long temporal horizons.

Figure 10. Momentum thickness of DNS, no-model and learned model simulations, evaluated based on the streamwise averages.

5. Spatially developing planar mixing layers

In contrast to the temporally developing mixing layers investigated in § 4, the spatially developing counterpart features a fixed view on a statistically steady flow field, which introduces a new set of challenges to the learning task. While the main difficulty in previous transient simulations was the modelling of an evolving range of turbulent scales, the statistically steady nature of the spatially developing mixing layer requires a reproduction of the turbulent statistics in its own statistically steady state. This in turn necessitates long-term accuracy and stability.

Spatially mixing layers develop from an instability in the shear layer. This instability is driven by a disturbance at the inlet, whose characteristics have great effect on the mixing layer growth (Ho & Huang Reference Ho and Huang1982). In a simulation environment, these disturbances are realised by a Dirichlet inlet boundary condition, where temporally varying perturbations are added to a steady mean flow profile. As proposed by Ko, Lucor & Sagaut (Reference Ko, Lucor and Sagaut2008), a suitable inlet condition including perturbations can be written as

(5.1)\begin{equation} u_{{in}}(y,t) = 1 + \frac{\Delta U}{2}\tanh(2y) + \sum_{d=1}^{N_d}\epsilon_d(1-\tanh^{2}(y/2))\cos(K_d y)\sin(\varOmega_dt), \end{equation}

where the number of perturbation modes $N_d=2$ holds for our simulations. Furthermore, we used inviscid wall conditions for the two $y$-normal spanwise boundaries, and the outflow boundary was realised by a simple Neumann condition with a stabilising upstream sponge layer. For all simulations, we set the characteristic velocity ratio $\Delta U=1$ and the vorticity thickness to $\delta _\omega =1$. The vorticity-thickness Reynolds number is set to $Re_{\delta _\omega }={\Delta U \delta _\omega }/{\nu } = 500$. To generate the DNS dataset, this set-up was discretised by a uniform grid with $(N_x,N_y)=(2048,512)$ resolving the domain of size $(L_x,L_y)=(256,64)$. The timesteps were chosen such that CFL$=0.3$ and the temporal evolution was run for $7$ periods of the slowest perturbation mode $i=2$ to reach a statistically steady state, before subsequent frames were entered into the dataset. A further $28$ periods of the slowest perturbation mode were simulated to generate $32\ 000$ samples of the statistically steady state. The training dataset consists of $5$ such simulations with different perturbations, as summarised in table 4. A downsampling ratio of $8\times$ in space and time was again chosen for the learning set-up. The input to the network was set to include only the main simulation frame without the sponge layer region.

Table 4. Perturbation details for the inlet condition of training and test datasets.

Our best performing model applied the turbulence loss $\mathcal {L}_{T}$, with the loss factors set to $(\lambda _2 , \ \lambda _E, \ \lambda _{S}, \ \lambda _{MS})= (50, \ 0.5, \ 2, \ 0.5)$, and an unrollment of $60$ solver steps. The timespan covered by these $60$ solver steps is comparable to a full period of the slowest perturbation mode. Using the roll-up frequency of the spatial mixing layer as the basis for the time scale $t_{f_\omega }=1/f_\omega$, $60$ solver steps unroll $0.85t_{f_\omega }$. As we detail in the following, our test metrics show that this approach of unrolling roughly one integral time scale yields the best results.

First, we evaluate the influence of unrollment in this test case. Once again, we show comparisons with additional set-ups, the parametric details of which can be found in table 5. Similar to the temporal mixing layer, the $60$ step model was trained using a gradient stopping technique. A $30$-step back-propagation subrange performed best again by maintaining long-term information while avoiding instabilities in the gradient calculation. This model is described as $\textrm {NN}_{60,\mathcal {L}_{T}}$ in this section. Details regarding the method are explained in § 6. The table shows that the simulation with the 60-step neural network outperforms the no-model baseline by an order of magnitude. For these evaluations, we assessed the model capabilities by running a CNN-corrected forward simulation. This simulation was initialised with a downsampled frame from the DNS test dataset in its fully developed state. This test dataset is generated with different inflow conditions, where the inlet forcing lies outside of the training range, making these evaluations an out-of-sample generalisation test. The variation in inlet forcing affects the location and intensity of the mixing layer roll-up and vortex merging. The simulation was run for $5000 \Delta t$, or $36$ periods of the slowest perturbation mode in order to obtain data from a statistically stable state. Despite this time frame being orders of magnitude longer than what is seen by the models at training time, the 60-step model retains a stable simulation that closely matches the behaviour of the DNS reference. Interestingly, this longer unrollment, of the order of one integral time scale, is crucial to arrive at a stable model. The models trained with shorter unrollment exhibit various degrees of spurious oscillations, especially the $10$-step model. These oscillations most likely originate from slight deviations in turbulent structures (e.g. vortex roll-up) inferred by the network. Since short unrollment models have never seen any further development of these self-exited structures, applying said models eventually causes even stronger unphysical oscillations downstream. As before, we omit purely data-driven models trained with pre-computed simulation states. These produce undesirable solutions within a few timesteps of simulating the test cases. The vorticity visualisations after half a period of the slowest perturbation mode ($70 \Delta t$) and after four periods or one flow through time ($600 \Delta t$) are shown in figures 11(a) and 11(b), and compared with DNS and the no-model simulation. The early evaluation in figure 11(a) reveals a severe loss of detail in the no-model simulation, even after a short time horizon. Over this timespan, the learned model achieves a close visual reproduction. Additionally, the later vorticity heat map in figure 11(b) shows a delayed roll-up in the no-model simulation, whereas the learned model maintains the roll-up location and shows improved accuracy. This behaviour is clarified by the Reynolds-averaged properties of the simulations, for which resolved Reynolds stresses and turbulence kinetic energies were calculated on the basis of the respective statistically steady simulations. As shown in figure 12, the no-model statistics severely deviate from the targeted DNS. In contrast, the corrective forcing inferred by the trained models approximates these statistics more accurately. The delayed roll-up of the no-model simulation and the improvement of the modelled ones is connected to the Reynolds stresses. The Reynolds stresses indicate turbulent diffusion of momentum, and figure 12 shows that the CNN learned to encourage turbulent fluctuations at the start of the mixing layer. The fluctuations trigger the shear layer instability and feed the roll-up, with decisive implications for the downstream development of the mixing layer. Especially the long unrollment of $60$ steps benefits the model performance. Evaluations at locations downstream of the initial roll-up see the accuracy of $10$ and $30$ step models deteriorate in direct comparison with the $60$-step model.

Figure 11. Vorticity heat maps of the spatial mixing layer simulations at (a) $t = 70\Delta t$, and (b) $t = 600\Delta t$, on the test dataset.

Figure 12. Comparison of downsampled DNS, no-model and learned model simulations with respect to Reynolds-averaged resolved turbulence kinetic energy (a) and Reynolds stresses (b).

Table 5. Model details for unrollment study; MSE with respect to DNS from test data at $t_e = 1000\Delta t$.

These observations regarding the Reynolds stresses extend to the resolved turbulence kinetic energies (figure 12), where the same turbulent fluctuations yield an accurate reproduction of the DNS. The learned models are not limited to a specific spatial scale, but precisely match the DNS on all turbulent scales when comparing the centreline kinetic energy spectra in figure 13.

Figure 13. Centreline kinetic energy spectra for downsampled DNS, no-model and learned model simulations.

The evaluations of vorticity and momentum thickness in figures 14(a) and 14(b) capture a delayed mixing layer development. Especially early stages of the mixing layer immediately after the first roll-up are modelled inaccurately. While all models show this behaviour, the delay in terms of momentum thickness is more pronounced for the long unrollment $60$-step model. On the contrary, the roll-up inaccuracy results in a noticeable offset in the vorticity thickness around $x/\delta _{\omega _0}=100$ for all models, but the $60$-step model performs best further downstream by recovering the DNS behaviour. This recovery is lacking in the $10$- and $30$-step models, causing the evaluation of Reynolds stresses at $x=192\Delta x$ (figure 12) to exhibit large discrepancies between DNS and learned model simulation for these models, with the notable exception of the $60$-step model. Note, however, that, despite not being capable of exactly reproducing the entire mixing layer up to the finest detail, the learned models still greatly outperform a no-model simulation. Momentum thickness evaluations show beneficial results for the models trained with shorter unrollments. Due to the definition of momentum thickness as an integral quantity over the shear direction, an increase in this quantity is caused by strong deviations from the initial step profile of the mixing layer. While the integral values for the momentum thickness of the $10$- and $30$-step models are close to the DNS, the underlying turbulence fluctuations causing these values are not accurate compared with the DNS, which can be seen in turbulence kinetic energy and Reynolds stress evaluations in figure 12. Considering these results jointly, we draw the conclusion that the $60$-step model yields the best performance.

Figure 14. Vorticity and momentum thickness of the downsampled DNS, no-model and learned model simulations.

Additionally, the evaluations show the benefits of training through multiple unrolled steps. The 10-step model develops instabilities after $500\Delta t$, which is equivalent to one flow-through time. From this time on, the learned model only sees self-exited instabilities in the mixing layer. This constitutes an extrapolation with respect to the temporal unrollment, as well as with respect to the inlet perturbation due to the use of a test dataset. This in turn can cause spurious oscillations and thus a deterioration of solution quality. The $30$-step model shows this behaviour to a lesser extent and generates a stable, statistically steady progression of the mixing layer for this case of temporal extrapolation. Even better behaviour is only achieved by the $60$-step model. It practically eliminates the instabilities seen in other models.

While previous evaluations showcased the stability improvements gained by training through multiple solver steps, another benefit of this approach relates to the temporal fluctuations in DNS training data. As visualised in figure 15, only some of the interactions between CNN and these temporal oscillations are covered in a training iteration. Consequently, the training loop imposes a high-pass cutoff on the observed frequencies that directly depends on the number of unrolled solver steps. To extract the temporal features that our models learned from the training dataset, we calculate the power spectral density of the velocity fields at sampling point $(x,y)=(160,0)$ on training data. The sampling timespan for the learned models starts after one flow-through time and stops after the next four flow-through times passed. The resulting power spectral densities are compared with a long-term evaluation of the DNS data, and a relative error between the spectra is computed. The results are shown in figure 15 and support the following observations. Firstly, all learned models can capture the discrete nature of the dominant frequencies quite well. In particular, the $60$-step model shows good approximation to the DNS evaluation. In contrast, the no-model does not match the DNS characteristics. Secondly, the relative error of the power spectra generated by the $60$-step model is substantially lower for all but the highest frequencies. Since the $30$- and $10$-step models only saw the interaction with fine scales during their training, these models perform worse at the lower frequencies, which results in higher relative errors for the relatively low vortex roll-up and vortex merging frequencies. These features operate at of the order of one integral time scale and are better resolved by 60 unrolled steps.

Figure 15. Power spectral density (PSD) of velocity fluctuations over time at sampling point $(x,y)=(192\Delta x,0)$ based on the training dataset for DNS, no-model and learned model simulations at top; bottom figure displays the relative error of the power densities over frequencies, accumulated for both velocity components; frequencies to the right of a dotted vertical line are fully enclosed in a training iteration; vertical lines correspond to (60, 30, 10) unrolled steps from left to right.

6. Gradient back-propagation

Our evaluations on temporally and spatially developing mixing layers show significant performance gain by longer unrollment times, with the best accuracy given by a $60$-step model. However, long unrollments can cause stability problems. Repeated applications of neural networks are known to be problematic during training, where exploding or diminishing gradients can significantly deteriorate the quality of gradients (Pascanu, Mikolov & Bengio Reference Pascanu, Mikolov and Bengio2013). To avoid this, we utilise a custom version of the gradient stopping technique: instead of discarding gradients generated by some (earlier) simulation steps, we split the gradient back-propagation into individually evaluated subranges. In other words, the training still exposes long temporal unrollments and preserves the gradient influence of all steps, but does not propagate gradients back to the first application of the network model. We use $60$-step models to study model accuracy with respect to the length of these back-propagation subranges on a range of $10$, $20$, $30$ and $60$ backward steps. We will use the notation $\textrm {NN}_{m\unicode{x2013} g}$ with two numbers $m$ and $g$, where $m$ describes the number of unrolled forward steps, and $g$ represents the length of the subranges for which gradients are calculated individually. In practice, this means that gradients of a $60$-$20$ model are only back-propagated through $3$ non-overlapping sections of $20$ steps each. The $60$-$60$ model recovers the standard differentiable training procedure used for previous models.

This procedure was applied to temporally and spatially developing mixing layers. Details of the trained models are found in tables 6 and 7. Note that the training of the $\textrm {NN}_{60\unicode{x2013} 60,\mathcal {L}_{T}}$ was not stable for the temporal mixing layer case, which we attribute to unstable gradients in the optimisation. In contrast, the subrange gradient models are stable during training. Additional evaluations of Reynolds stresses and turbulence kinetic energy for the temporal mixing layer indicate no performance differences between these models, as shown in figure 16. We thus conclude that the method of subrange back-propagation makes the training of the 60-step model possible, but also that the model performance on the temporal mixing layer was already saturated by the $30$-step model, as previously mentioned in § 4. The $\mathrm {NN}_{60\unicode{x2013} 30,\mathcal {L}_{T}}$ model was used in the evaluation in § 4.

Figure 16. Comparison of DNS, no-model and $60$-step model simulations with respect to resolved turbulence kinetic energy (a), and Reynolds stresses (b).

Table 6. Temporal mixing layer; model details for $60$-$X$ models; MSE with respect to DNS from test data at $t_e = 512\Delta t$; training of $\textrm {NN}_{60\unicode{x2013} 60,\mathcal {L}_{T}}$ is unstable.

Table 7. Spatial mixing layer; model details for $60$-$X$ models; MSE with respect to DNS from test data at $t_e = 1000\Delta t$.

The spatial mixing layer models are evaluated on vorticity snapshots in figure 17, turbulence kinetic energy and Reynolds stresses in figure 18, as well as vorticity and momentum thickness in figures 19. These results indicate that there is a optimal number of consecutive back-propagation steps of around $20$ to $30$, where the optimisation gradients contain long-term information while still maintaining a good quality that is unaffected by risks of recurrent evaluation. The $\mathrm {NN}_{60\unicode{x2013} 20,\mathcal {L}_{T}}$ and $\mathrm {NN}_{60\unicode{x2013} 30,\mathcal {L}_{T}}$ models achieve best performance on all metrics except for the momentum thickness. We attribute the larger values of momentum thickness to some spurious oscillations exhibited by the $\mathrm {NN}_{60\unicode{x2013} 10,\mathcal {L}_{T}}$ and $\mathrm {NN}_{60\unicode{x2013} 60,\mathcal {L}_{T}}$ models. The $\mathrm {NN}_{60\unicode{x2013} 30,\mathcal {L}_{T}}$ model was used in earlier unrollment evaluations in § 5.

Figure 17. Vorticity comparison of $60$-step models on spatial mixing layer simulations at $t=700\Delta t$ on the test dataset.

Figure 18. Comparison of downsampled DNS, no-model and learned model simulations with respect to Reynolds-averaged resolved turbulence kinetic energy (a) and Reynolds stresses (b).

Figure 19. Vorticity and momentum thickness of the downsampled DNS, no-model and learned model simulations.

Another potential problem could be caused by training towards matching frames separated by long timespans. Turbulent flows could potentially loose correlation to the reference data over long temporal horizons, which would render this learning approach driven by simulated DNS data inapplicable. The unrollment times in this paper are, however, far from reaching an uncorrelated state. As shown in the previous evaluations, the $60$-step models perform better than their $30$-step counterparts, indicating that there is additional information provided by unrolling $60$ steps. This shows that the unrolled temporal horizons are far from exhibiting flow decorrelation. Further experiments with even longer unrollments on the spatial mixing layer revealed that no improvement is achieved beyond $60$ steps in this case. Figure 20 depicts selected evaluations of a 120-step model, which lack improvements over the $60$-step counterpart. While the $120$-step model gains accuracy in early upstream cross-sections, the mixing layer shift downstream of the first roll-up is worse in direct comparison.

Figure 20. Selected evaluations of the 120-step model with vorticity snapshots in (a), Reynolds stresses in (b) and vorticity thickness in (c).

We also investigated yet longer horizons (180 and 240 steps), but these runs saw a reduced accuracy with respect to some of the evaluations. One explanation is that the flow field is uncorrelated to the DNS data for these long horizons, leading to a diffused learning signal. If the loss was computed on late, uncorrelated frames, we would expect generated gradients to resemble random noise. While earlier frames would still provide valuable information, the random noise from these later frames could prevent the learning of precise corrections. In addition, the longer runs used the same set of hyperparameters as determined for the shorter unrollments, the long horizon runs could also profit from a broader hyperparameter search.

In this section, we have identified gradient instabilities as the main problem when unrolling long temporal horizons. We have introduced a gradient splitting technique that stabilised the training procedure. This is done by splitting the gradient calculation into non-overlapping subranges. For the studied set-ups and $60$-step models, a split into two subranges of $30$ steps each performed best. One can conclude that longer unrollments pay off in terms of modelling accuracy up to a certain saturation point. In our simulations this saturation point lies at approximately $60$ steps, which coincides with the integral time scales of the respective scenarios. Unrolling beyond that saturation point is possible, but leads to increased computational effort and may require special treatment, such as a further tuning of the hyperparameters.

7. Computational performance

The development of turbulence models is ultimately motivated by a reduced computational cost, which facilitates numerical simulations in flow scenarios where a DNS is prohibitively expensive. Preceding sections have outlined the corrective capabilities of our learned models. We now seek to put these improvements into perspective by studying the computational cost of our learned models at the inference time. For all of our performance evaluations, an Intel Xeon E5-1650 CPU and a Nvidia GTX 1080Ti GPU are used. We use the computational set-ups from our model evaluation runs on test data in the isotropic turbulence, TML and SML cases in §§ 3, 4 and 5, respectively.

Exactly as before, an $8\times$ scaling factor is deployed on both the spatial resolution and timestep size. We then run the simulations until the time $t_e=1000\Delta t$ is reached, while recording the required computational time for each timestep. The results are summarised in table 8, where the total simulation time as well as per-timestep values are listed. We also assess the computational cost of a no-model simulation that matches the performance of our models.

Table 8. Computational performance comparison over $t_e=1000\Delta t$ for the used flow scenarios, isotropic decaying turbulence (IDT), temporal mixing layer (TML) and spatial mixing layer (SML); MSE values are evaluated on the velocity field at $500\Delta t$; training time on one GPU.

The resulting data show that the neural network incurs only a negligible cost of approximately $10\,\%$ in comparison with no-model simulations at the same resolution. The learned models clearly outperform the no-model variants in terms of MSEs, and incur only a fraction of the computational cost required for the DNS variants.

In addition, we provide the temporal evolution of the MSE evaluated on resolved turbulence kinetic energies for all three scenarios in figure 21. From this evaluation, we conclude that our method consistently outperforms simulations with a $2\times$ higher resolution in spatial and temporal dimensions. Additionally, we found our learned models to often be on-par with $4\times$ higher resolved simulations, e.g. in the first half of the temporal mixing layer case. On the basis of the clock times from table 8, this corresponds to a speed up of $3.3$ over $2\times$ isotropic turbulence simulations. For the mixing layer cases, the hybrid model on average resembles the performance of $3\times$ reference simulations, which corresponds to a speed-up of $7.0$ for the temporal and $3.7$ for the spatial mixing layer. For the former, our model even closely matches the performance of a $4\times$ simulation for several hundred timesteps, which represents a speed up of $14.4$.

Figure 21. Similarity evolutions over time measured by the MSE on resolved turbulence kinetic energy for randomised turbulence simulations (a), TML simulations (b) and SML simulations (c).

While other works have reported even larger performance improvements (Kochkov et al. Reference Kochkov, Smith, Alieva, Wang, Brenner and Hoyer2021), we believe that our measurements are representative of real-world scenarios with higher-order solvers. Asymptotically, we also expect even larger payoffs for the high-resolution, three-dimensional simulations that are prevalent in real-world applications.

Naturally, the training of each neural network requires a substantial one-time cost. In our case, the network took 3 to 10 days of training, depending on the individual problem set-up. The required GPU hours for the best-performing models are listed in table 8. The longer unrolled temporal horizons and larger domain increase the required training time for the spatial mixing layer. For the three used set-ups, these training times are equivalent to [120, 118, 22] DNS solves of full length as used in the dataset calculation. However, under the assumption that the learned turbulence model can be employed by multiple users in a larger number of simulations to produce new outputs, this cost will quickly amortise. Especially the most complex spatial mixing layer case shows a favourable relation of training cost to simulation speed up. Additionally, a successful application of this approach to three-dimensional turbulence would make training cheaper in relation to DNS and speed ups larger, due to the scaling through an additional spatial dimension. It is worth noting that our comparisons are based on GPU solvers, and performance is likely to vary on CPU or mixed solvers, where parts of the computation are CPU-based and communication overheads could deteriorate gains.