2 Governing equations

We proceed by outlining our Kraichnan turbulence test case which is a simplified prototype for geophysical flow phenomena (Pearson et al. Reference Pearson, Fox-Kemper, Bachman and Bryan2017). The governing equations of motion for Kraichnan turbulence are given by the two-dimensional Navier–Stokes equations in a periodic domain. The non-dimensionalized version of these equations may be expressed in the vorticity ( $\unicode[STIX]{x1D714}$ ) and streamfunction ( $\unicode[STIX]{x1D713}$ ) formulation as (San & Staples Reference San and Staples2012),

(2.1) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}t}}+J(\unicode[STIX]{x1D714},\unicode[STIX]{x1D713})={\displaystyle \frac{1}{Re}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D714},\\ x,y\in [0,2\unicode[STIX]{x03C0}],\quad t\in [0,4],\end{array}\right\}\end{eqnarray}$$

where we define the Jacobian (or the nonlinear term as)

(2.2) $$\begin{eqnarray}J(\unicode[STIX]{x1D714},\unicode[STIX]{x1D713})={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}x}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}y}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}{\unicode[STIX]{x2202}y}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}x}},\end{eqnarray}$$

and the conservation of mass is enforced by

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}=-\unicode[STIX]{x1D714}.\end{eqnarray}$$

A measure of multi-scale behaviour in this system is given by the Reynolds number ( $Re$ ). A high value of $Re$ combined with a coarse-grid projection of these equations results in insufficient support for the finest structures in the flow evolution, leading to noise accumulation at grid cutoff and potential floating point overflow of the numerical evolution of this problem. A sufficiently coarse-grained representation of the governing equations introduced previously are given by the LES governing equations

(2.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}{\displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D714}}}{\unicode[STIX]{x2202}t}}+J(\bar{\unicode[STIX]{x1D714}},\bar{\unicode[STIX]{x1D713}})={\displaystyle \frac{1}{Re}}\unicode[STIX]{x1D6FB}^{2}\bar{\unicode[STIX]{x1D714}}+\unicode[STIX]{x1D6F1},\\ \unicode[STIX]{x1D6FB}^{2}\bar{\unicode[STIX]{x1D713}}=-\bar{\unicode[STIX]{x1D714}},\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D6F1}$ may be assumed to be the perfect closure given by

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}=J(\bar{\unicode[STIX]{x1D714}},\bar{\unicode[STIX]{x1D713}})-\overline{J(\unicode[STIX]{x1D714},\unicode[STIX]{x1D713})}.\end{eqnarray}$$

When adequately simulated, the decaying Kraichnan turbulence test cases result in the classical $k^{-3}$ scaling of the energy spectra (Kraichnan Reference Kraichnan1967). We focus on two competing ideologies for approximating closure. The first is given by the functional hypothesis and may be expressed as

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}=\unicode[STIX]{x1D708}_{e}\unicode[STIX]{x1D6FB}^{2}\bar{\unicode[STIX]{x1D714}},\end{eqnarray}$$

where the Smagorinsky approximation to the eddy viscosity $\unicode[STIX]{x1D708}_{e}$ is given by

(2.7) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D708}_{e}=(C_{s}\unicode[STIX]{x1D6FF})^{2}|\bar{S}|,\\ |\bar{S}|=\sqrt{4\Big({\displaystyle \frac{\unicode[STIX]{x2202}^{2}\bar{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y}}\Big)^{2}+\Big({\displaystyle \frac{\unicode[STIX]{x2202}^{2}\bar{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x2202}x^{2}}}-{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\bar{\unicode[STIX]{x1D713}}}{\unicode[STIX]{x2202}y^{2}}}\Big)^{2}},\end{array}\right\}\end{eqnarray}$$

where it is very common to consider the filter length scale $\unicode[STIX]{x1D6FF}$ as the representative mesh size. The static specification of $C_{s}$ might yield over-dissipative or under-dissipative results (e.g. see Canuto & Cheng Reference Canuto and Cheng1997; Vorobev & Zikanov Reference Vorobev and Zikanov2008; Cushman-Roisin & Beckers Reference Cushman-Roisin and Beckers2011). A successful application of this closure necessitates a dynamic calculation of the Smagorinsky coefficient $C_{s}$ that requires the specification of a test filter and a spatial averaging for stabilized deployment. This approach is the well-known dynamic Smagorinsky (DS) closure (Germano et al. Reference Germano, Piomelli, Moin and Cabot1991; Lilly Reference Lilly1992) and its two-dimensional abstraction for Kraichnan turbulence has been presented by San (Reference San2014) by showing that $C_{s}$ dynamically varies between 0.1 and 0.2 in agreement with previously reported values in LES literature (e.g. see also Maulik & San (Reference Maulik and San2017a ) for a modular approach in various functional forms).

A competing ideology is given by the structural (or scale-similarity) hypothesis which assumes that the LES equations are projections of the Navier–Stokes equations to a smoother space where an inverse-filtering operator may be utilized to recover the finer scales that are lost. Mathematically,

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}=J(\bar{\unicode[STIX]{x1D714}},\bar{\unicode[STIX]{x1D713}})-\widetilde{J(\unicode[STIX]{x1D714}^{\ast },\unicode[STIX]{x1D713}^{\ast })},\end{eqnarray}$$

where $\unicode[STIX]{x1D714}^{\ast }$ and $\unicode[STIX]{x1D713}^{\ast }$ are approximately deconvolved variables obtained through an inverse-filtering procedure and a Gaussian-type filter kernel (given by $G(a)=\tilde{a}$ with $G$ being the filter kernel applied on a field variable $a$ and $\tilde{a}$ being its low-pass spatially filtered counterpart) is an approximation of the projection to the LES space. However, these techniques are limited due to the underlying assumption of isomorphism between the LES and the Navier–Stokes equations (Germano Reference Germano2015). In practice, this implies that structural hypotheses are appropriate only if finer structures are sufficiently well resolved on a particular grid. As such, this diminishes their benefit for practical flows where grid cutoff wavenumbers are generally much smaller than the largest wavenumber in the flow. The breakdown of structural closures manifests itself in the form of stability issues. For this reason, many successful closure deployments utilize linear combinations of structural and functional models (Habisreutinger et al. Reference Habisreutinger, Bouffanais, Leriche and Deville2007). In this work, we implement approximate deconvolution (AD) (Stolz & Adams Reference Stolz and Adams1999; San et al. Reference San, Staples, Wang and Iliescu2011) which utilizes an iterative application of the trapezoidal filter kernel for inversion of filtered grid quantities, and utilize three iterative resubstitutions to deconvolve our grid-resolved variables.

3 Machine learning

We now turn to the procedure of utilizing DNS data for learning when to switch between one of three closure scenarios. Of these three options, two are given by the choice of the functional hypothesis and AD. The third option is that of a no-model scenario where our learning determines that closure modelling is unnecessary. The third scenario is retained since there is a possibility of localized areas in a flow having adequate grid support so that the contributions of the sub-grid scale become negligible. This switching between scenarios is spatio-temporally dynamic. Before proceeding, we note that the functional deployment eschews the dynamic procedure and simply sets an arbitrarily large value of $C_{s}=1.0$ for the calculation of the eddy-viscosity kernel given by (2.7) as utilized in our training. We proceed by outlining our learning strategy through the utilization of DNS data. Five snapshots of DNS data at $Re=32\,000$ and at $2048^{2}$ degrees of freedom (from 40 000 available snapshots) are utilized to compute the grid-filtered variables (herein denoted by FDNS) at $256^{2}$ degrees of freedom through the application of a spectral cutoff filter (i.e. only retaining the wavenumbers with respect to the coarse-grained LES grid support). Perfect closure values ( $\unicode[STIX]{x1D6F1}$ ) are then obtained (the reader is directed to Maulik et al. (Reference Maulik, San, Rasheed and Vedula2019) for details related to the calculation of these quantities). Figure 1 visually quantifies the effect of the spectral domain filtering where the FDNS of a snapshot of vorticity is shown. We then introduce the a priori eddy viscosity given by

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D708}_{e}^{a}={\displaystyle \frac{\unicode[STIX]{x1D6F1}}{\unicode[STIX]{x1D6FB}^{2}\bar{\unicode[STIX]{x1D714}}}},\end{eqnarray}$$

where all the terms on the right-hand side of the above equation are available through calculation from the DNS snapshots. The a priori eddy viscosity is centred at a value of zero (corresponding to a region where closure modelling is unnecessary) and has tails in the negative and positive directions. A core component of the hypothesis in this work stems from the fact that structural hypotheses are not limited to positive eddy-viscosity predictions alone. The reader may note that models utilizing the functional hypothesis always lead to positive eddy viscosities. The a priori eddy viscosities calculated from the DNS data are then projected onto a Gaussian distribution where values lying within a distance of 1 % of the standard deviation from the mean (of zero) are labelled as those requiring no closure (due to the low strength of the a priori eddy viscosity). Values lying beyond this range are labelled as functional or structural, depending on whether they are positive or negative, respectively. This information is encoded in one-hot labelling for a classification deployment and a corresponding schematic for this hypothesis segregation is shown in figure 2. It is observed that a large portion of the available data lies within the first standard deviation of the mean eddy viscosity. This leads to the potential of turbulence modelling classification being considered from an outlier detection point of view. A factor which motivates the choice of the Gaussian distribution is the nature of the decaying Kraichnan turbulence. However, we note that machine learning algorithms are also capable of classifying data belonging to complex distributions and that this hypothesis segregation may be tuned for better accuracy. Also, the choice of the 1 % hyper-parameter is also motivated by observing a posteriori training accuracies where it is noticed that a relatively simple architecture (mentioned next) is efficiently able to discern the varying hypothesis. Values greater than 1 % for model delineation led to reduced learning accuracies indicating that a physical delineation potentially exists in this projection and categorization. Essentially, categorizing closure bins for values greater than 1 % leads to poorer training and validation classification accuracies during the machine learning phase. This could indicate that a clear distinction between modelling requirements is lost through this segregation. Further study for adding complexity to the hypothesis segregation is thus, a necessity.

Figure 1.

Visualization of the effect of Fourier cutoff filtering with DNS at $N^{2}=2048^{2}$ (a) and corresponding FDNS at $N^{2}=256^{2}$ (b).

Each label for the a priori eddy viscosity is also associated with an input kernel of grid-resolved quantities. This kernel is given by a stencil of 9 inputs each of vorticity and streamfunction (for a total of 18 input variables). These 9 inputs of each field are given by a query of the field quantity at a point on the coarse grid, the 4 adjacent points in each dimension and the 4 diagonally adjacent points. Each sample of our training data thus consists of 18 inputs of vorticity and streamfunction and outputs given by one-hot labels for the choice of closure modelling strategy. In this article, we have leveraged the fact that the means of vorticity and streamfunction are both very close to zero and do not necessitate normalization. In addition, the non-dimensionalized formulation of the governing equations implies that our inputs are all dimensionless. However, we note that for practical deployments of any local-kernel-based machine learning queries, grid-resolved quantities must be normalized and non-dimensionalized.

Figure 2.

Data segregation for one-hot labelling. The a priori eddy viscosities are projected onto a Gaussian distribution where data beyond 1.0 % of the standard deviation are labelled as requiring structural (if negative) or functional (if positive) modelling. The remaining data points are classified as no-model cases.

Mathematically, if our input vector $\boldsymbol{p}$ resides in a $P$ -dimensional space and our desired output $\boldsymbol{q}$ resides in a $Q$ -dimensional space, this framework establishes a map $\mathbb{M}$ as follows:

(3.2) $$\begin{eqnarray}\mathbb{M}:\{p_{1},p_{2},\ldots ,p_{P}\}\in \mathbb{R}^{P}\rightarrow \{q_{1},q_{2},\ldots ,q_{Q}\}\in \mathbb{R}^{Q}.\end{eqnarray}$$

Accordingly, the framework utilized in this article leads to the following relation:

(3.3) $$\begin{eqnarray}\mathbb{M}:\{\boldsymbol{p}\}\in \mathbb{R}^{18}\rightarrow \{P(\boldsymbol{q}|\boldsymbol{p})\}\in \mathbb{R}^{3},\end{eqnarray}$$

where our input and output spaces are given by

(3.4) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{p}_{i,j}=\{\bar{\unicode[STIX]{x1D714}}_{i,j},\bar{\unicode[STIX]{x1D714}}_{i,j+1},\bar{\unicode[STIX]{x1D714}}_{i,j-1},\ldots ,\bar{\unicode[STIX]{x1D714}}_{i-1,j-1},\bar{\unicode[STIX]{x1D713}}_{i,j},\bar{\unicode[STIX]{x1D713}}_{i,j+1},\bar{\unicode[STIX]{x1D713}}_{i,j-1},\ldots ,\bar{\unicode[STIX]{x1D713}}_{i-1,j-1}\},\\ \boldsymbol{q}_{i,j}=\{P(\unicode[STIX]{x1D6F1}_{i,j}^{k}|\boldsymbol{p}_{i,j})\},\end{array}\right\}\end{eqnarray}$$

where $i,j$ refer to the spatial indices on the coarse grid (i.e. the point of deployment) and $k$ refers to the choice of closure scenario (i.e. structural, functional or no closure). We note here that the choice of the local stencil for ANN query reflects the discretization of the governing equations (with second-order accurate stencils requiring a $\pm 1$ query) and the use of the trapezoidal filter in AD. Also, note that our choice of input space is given by raw variable queries rather than derivatives (or other such engineered terms). This is motivated by an aversion to specify bias towards any particular quantity that may otherwise by learned implicitly by the network. However, we note that the classification workflow may benefit significantly from the inclusion of a feature engineering step prior to optimization. This is a subject of ongoing investigation.

Our optimal map $\mathbb{M}$ is then trained by the following loss function

(3.5) $$\begin{eqnarray}E(\boldsymbol{w})=-\mathop{\sum }_{n=1}^{N}\mathop{\sum }_{k=1}^{K}\{t_{nk}\log (y_{nk})+(1-t_{nk})\log (1-y_{nk})\},\end{eqnarray}$$

where $\boldsymbol{w}$ are the tuneable weights and biases of the network, $N$ is the total number of samples and $K=3$ is the total number of closure scenarios. Here, $t_{nk}$ refers to the target (or true) label of class $k$ and sample $n$ and $y_{nk}$ refers to its corresponding prediction. Note that one-hot encoding ensures that $t_{nk}$ values are always binary (Bishop Reference Bishop2006). For reference, our architecture is trained using the open-source deep learning software Tensorflow and is optimized with the use of ADAM, a popular gradient-descent-based optimizer.

Our learning architecture is given by a 5 hidden-layer deep neural network with 40 neurons each for calculating the conditional probabilities of the three closure scenarios pointwise in space and time. The hidden-layer neurons employ a rectified-linear activation and the output layer gives us softmax probabilities for the three classes. The scenario with the highest conditional probability is then deployed for model classification i.e.

(3.6) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{i,j}^{ML}=\unicode[STIX]{x1D6F1}_{i,j}^{k}\text{ s.t. }\text{argmax}_{k}P(\unicode[STIX]{x1D6F1}_{i,j}^{k}|\boldsymbol{p}_{i,j}),\end{eqnarray}$$

where $\unicode[STIX]{x1D6F1}_{i,j}^{ML}$ refers to the machine-learning-based turbulence model computation at a point. In the case of model blending, the conditional probabilities for closure scenarios are used to find a linear combination of the standard Smagorinsky and the AD closures. In other words,

(3.7) $$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{i,j}^{ML}=P(\unicode[STIX]{x1D6F1}_{i,j}^{AD}|\boldsymbol{p}_{i,j})\unicode[STIX]{x1D6F1}_{i,j}^{AD}+P(\unicode[STIX]{x1D6F1}_{i,j}^{SM}|\boldsymbol{p}_{i,j})\unicode[STIX]{x1D6F1}_{i,j}^{SM},\end{eqnarray}$$

where $\unicode[STIX]{x1D6F1}_{i,j}^{AD}$ and $\unicode[STIX]{x1D6F1}_{i,j}^{SM}$ are AD and Smagorinsky predictions for the turbulence model at a point. We note that the same learning framework is deployed in these two conceptually different scenarios.

The framework is trained using the previously mentioned categorical cross-entropy error minimization for the one-hot encoded targets. A threefold cross-validation is utilized with a grid search for the number of layers (between 1 and 8) and number of neurons (between 10 and 100 at intervals of 10) to arrive at the optimal architecture mentioned previously. This optimal network is then trained for 2000 epochs to arrive at a classification accuracy of 79 % for training and approximately 68 % accuracy for validation. Convergence in validation loss was observed at around 1500 epochs as shown in figure 3. We note here that our validation data (amounting to one third of the total training data set) were not exposed to the network during gradient calculation in the back-propagation-based training procedure. Effectively, our learning is derived from two thirds of the total training data while our best model is chosen from that with the lowest validation loss. This is to ensure that the chances of network extrapolation are minimized. Multiple random initializations of our network weights led to similar training and validation behaviour as also seen in figure 3. Note one initialization is seen to take a longer time to converge as against the rest. The optimal learning is then deployed into a posteriori evolution of the Kraichnan turbulence test case where a pointwise closure deployment is performed for a variety of test cases. We also note that our labelled data are pre-processed to ensure that an equal number of samples are available from each classification regime to prevent our learning from prioritizing one outcome over the other two.

Figure 3.

Neural network training and validation loss for the proposed learning framework for one random weight initialization showing convergence at around 1500 epochs (a). The training (in black) and validation (in red) saturation behaviour for different initializations also indicated convergence to similar accuracies as shown on the right. The best model was chosen according to lowest validation loss for reduced over-fitting in forward deployments.

5 Conclusions and significance

In this article we have proposed a novel data-driven strategy to dynamically assess the utility of a turbulence modelling hypothesis in an LES framework. This strategy is built on the hypothesis that DNS data may be utilized to assess areas where structural or functional models may be more appropriate in an LES deployment. Our hypothesis segregation and subsequent training culminates in a learning that may deployed as a classifier of turbulence models at each point on the LES grid as well as a blending technique for balancing turbulence models with different dissipative strengths. When deployed as a classifier, our proposed framework may also predict a ‘no-model’ situation wherein no sub-grid source-term is deployed. When deployed as a blending mechanism, the learning linearly combines the AD and static Smagorinsky hypothesis premultiplied by their respective conditional probabilities to obtain another hybrid dissipation mechanism. Both frameworks utilize the same learning and are assessed through similar experiments a posteriori. Furthermore, a similar architecture can be constructed to provide hybrid intelligent implicit LES approaches by toggling between non-dissipative and dissipative numerical schemes, a topic we intend to investigate further in a future study.

We have rigorously assessed the deployment of our machine learning strategy through the utilization of a Kraichnan turbulence test case. Our assessments are made for Reynolds number values both within and outside that utilized in training to ensure that a generalizable turbulence closure has been developed. In addition, we have also assessed if the proposed closure can be deployed on a coarser grid than one it was trained for. The dissipative and scale-content capture of the proposed framework is compared to the AD and DS techniques through the use of kinetic-energy spectra, vorticity structure functions and time histories of $TKE$ and vorticity variance showing a dynamic dissipation akin to the DS. In particular, the statistical fidelity of the data-driven frameworks is seen to be inferior to the AD technique, which provides better estimates of the kinetic-energy spectra at lower wavenumbers and also provides most accurate estimates of the vorticity structure function. However, the focus on high-wavenumber noise attenuation leads to no grid cutoff error accumulation and the statistical results of the ML models are very similar to DS in all assessments. Also, it is observed that the data-driven closure (whether deployed as a classifier or a blending instrument) adequately captures the $k^{-3}$ scaling expected for the kinetic-energy spectra for the Kraichnan turbulence case and attempts to strike an optimal balance between the dissipative functional kernel and the noise-prone structural kernel. This behaviour is interesting as the model classifies solely between AD and the static Smagorinsky hypothesis indicating the extreme dissipation of the latter at $C_{s}=1.0$ is effectively alleviated by the spatio-temporal blending. Our closure, thus, attempts to blend the strengths of both modelling strategies to overcome their individual weaknesses while attempting to preserve trends from DNS.

In terms of future opportunities for this idea, the data-driven element of closure identification lends to the potential development of closures that may discern the physical characteristics of different flow scenarios. However, some challenges associated with progress in this research include considerations of invariance properties, which we have identified as a next step. In order to generalize adequately, Galilean invariance must be built in to the data-driven closure described in this study. This is a limitation of the current approach which does not explicitly enforce these considerations. For instance, research is underway to investigate the response of this system to invariant inputs of localized stencils (rather than raw field variables) to ensure that predictions are invariant as well. We may also look at augmenting our input data set through the use of transformations which make the classifier resistant to extrapolation (Ling et al. Reference Ling, Jones and Templeton2016a ). We note that the model blending idea proposed here can be devised to ensure invariance properties by choosing a suitable projection basis (Ling et al. Reference Ling, Kurzawski and Templeton2016b ). Our frameworks may then be assessed for transformed validation data. We expect an optimal formulation that relies on invariant inputs to map to a basis of invariant outputs for predicting an appropriate closure as a culmination of this study. The effect of input spaces (invariant or otherwise) on generalizability and accuracy must also be studied. This is an extension to the current formulation which takes a rudimentary approach of exposing the framework to all available grid-resolved information in the local stencil. Another major limitation of our study is the finite number of assessments in decaying two-dimensional turbulence where the vortex stretching mechanism is absent.

While the computational costs of the proposed framework have not been studied in detail, an efficient deployment of the proposed framework would need graphical processing unit integration of any practical CFD simulation. The latter would lead to efficient learning queries since all the spatial domain information would be available to the common memory. Overall, without performing any additional effort by optimizing our deployment code, the proposed ML approach adds computational overhead that is approximately 2.5 times greater than that of required for simulating the problem with the DS model due to arithmetic computations within 5 hidden-layer network architecture at each time step. In practice, the trade-off between the training accuracy (hence the number of layers and the number of neurons at each layer) and computational efficiency of the ML-based models should be carefully analysed. A priori hyper-parameter search might yield a deep network architecture for an accurate training representation, but it might provide a costly deployment model for forward simulations. This would suggest to perform feature selection and simplify network architecture as recommended by Wang et al. (Reference Wang, Luo, Li, Tan and Fan2018). Another future direction identified in this research is the exposure of different two-dimensional turbulence physics to the classification framework to identify if closure choices can also be influenced by the training data regime. Success in that regard would allow for ‘train and forget’ closures in problems that have unsteady physics that span fundamentally different turbulence modelling requirements.