A priori analysis on deep learning of subgrid-scale parameterizations for Kraichnan turbulence

Abstract

In the present study, we investigate different data-driven parameterizations for large eddy simulation of two-dimensional turbulence in the a priori settings. These models utilize resolved flow field variables on the coarser grid to estimate the subgrid-scale stresses. We use data-driven closure models based on localized learning that employs a multilayer feedforward artificial neural network with point-to-point mapping and neighboring stencil data mapping, and convolutional neural network fed by data snapshots of the whole domain. The performance of these data-driven closure models is measured through a probability density function and is compared with the dynamic Smagorinsky model (DSM). The quantitative performance is evaluated using the cross-correlation coefficient between the true and predicted stresses. We analyze different frameworks in terms of the amount of training data, selection of input and output features, their characteristics in modeling with accuracy, and training and deployment computational time. We also demonstrate computational gain that can be achieved using the intelligent eddy viscosity model that learns eddy viscosity computed by the DSM instead of subgrid-scale stresses. We detail the hyperparameters optimization of these models using the grid search algorithm.

Appendices

Appendix A: Derivation of the Smagorinsky model in 2D turbulence

From Eq. 5, the subgrid-scale stresses in 2D field can be written as

$$\begin{aligned} \tau _{ij}&= \overline{u_i u_j} - {\bar{u}}_i {\bar{u}}_j, \end{aligned}$$

(33)

$$\begin{aligned}&= \underbrace{\frac{1}{2}\tau _{kk}\delta _{ij}}_{k_{\mathrm{SGS}}\delta _{ij}} + \bigg (\underbrace{ \tau _{ij} - \frac{1}{2}\tau _{kk}\delta _{ij}}_{\tau _{ij}^d} \bigg ). \end{aligned}$$

(34)

The SGS stresses can be written as

$$\begin{aligned} \tau = k_{\mathrm{SGS}}I + \tau ^d, \end{aligned}$$

(35)

where $k_{\mathrm{SGS}}=\frac{1}{2}\tau _{kk}$ is called subgrid-scale kinetic energy (i.e., using the conventional summation notation with repeating indices, for example, $\tau _{kk} = \tau _{11} + \tau _{22}$, in 2D). In Smagorinsky model, we model the deviatoric (traceless) part of SGS stresses as

$$\begin{aligned} \tau _{ij}^d = -2\nu _\mathrm{e}{\bar{S}}_{ij}^d, \end{aligned}$$

(36)

where $\nu _\mathrm{e}$ is the SGS eddy viscosity, and ${\bar{S}}_{ij}$ is called resolved strain rate tensor given by

$$\begin{aligned} {\bar{S}}_{ij} = \frac{1}{2}\bigg ( \frac{\partial {\bar{u}}_i}{\partial x_j} + \frac{\partial {\bar{u}}_j}{\partial x_i} \bigg ), \end{aligned}$$

(37)

where we can write explicitly as follows

$$\begin{aligned} {\bar{S}} = \begin{bmatrix} \frac{\partial {\bar{u}}}{\partial x} &{}\quad \frac{1}{2}\bigg ( \frac{\partial {\bar{u}}}{\partial y} + \frac{\partial {\bar{v}}}{\partial x} \bigg ) \\ \frac{1}{2}\bigg ( \frac{\partial {\bar{v}}}{\partial x} + \frac{\partial {\bar{u}}}{\partial y} \bigg ) &{}\quad \frac{\partial {\bar{v}}}{\partial y} \end{bmatrix} . \end{aligned}$$

(38)

The trace of the ${\bar{S}}$ is zero owing to the continuity equation for incompressible flows. Therefore, ${\bar{S}}_{ij}^d = {\bar{S}}_{ij}$ and the Smagorinsky model becomes

$$\begin{aligned} \tau _{ij}^d = -2\nu _\mathrm{e}{\bar{S}}_{ij}. \end{aligned}$$

(39)

The eddy viscosity approximation computes $\nu _\mathrm{e}$ using the following relation

$$\begin{aligned} \nu _\mathrm{e} = C_k \varDelta \sqrt{k_{\mathrm{SGS}}}, \end{aligned}$$

(40)

where the proportionality constant is often set to $C_k = 0.094$, and $\varDelta$ is the length scale (usually grid size). The SGS kinetic energy $k_{\mathrm{SGS}}$ is computed with the local equilibrium assumption of the balance between subgrid-scale energy production and dissipation

$$\begin{aligned} {\bar{S}}:\tau + C_{\epsilon }\frac{k_{\mathrm{SGS}}^{1.5}}{\varDelta } = 0, \end{aligned}$$

(41)

where the first term in the above equation is dissipation flux, second term is production flux, and the production constant is often set to $C_{\epsilon }=1.048$. The double inner product operation  :  is given by

$$\begin{aligned} {\bar{S}}:\tau = {\bar{S}}_{ij}\tau _{ij}={\bar{S}}_{11}\tau _{11} + {\bar{S}}_{12}\tau _{12} + {\bar{S}}_{21}\tau _{21} + {\bar{S}}_{22}\tau _{22}. \end{aligned}$$

(42)

Substituting Eqs. 35 and 39 into Eq. 41, we get

$$\begin{aligned} {\bar{S}}:(k_{\mathrm{SGS}}I - 2C_k \varDelta \sqrt{k_{\mathrm{SGS}}}{\bar{S}}) + C_{\epsilon }\frac{k_{\mathrm{SGS}}^{1.5}}{\varDelta }&= 0, \end{aligned}$$

(43)

$$\begin{aligned} \sqrt{k_{\mathrm{SGS}}}\bigg ( \frac{C_{\epsilon }}{\varDelta }k_{\mathrm{SGS}} + \sqrt{k_{\mathrm{SGS}}} \underbrace{{\bar{S}}:I}_{{\bar{S}}_{ij} \delta _{ij} = 0} - 2C_k \varDelta {\bar{S}}:{\bar{S}} \bigg )&= 0, \end{aligned}$$

(44)

$$\begin{aligned} \frac{C_{\epsilon }}{\varDelta }k_{\mathrm{SGS}} - 2C_k \varDelta {\bar{S}}:{\bar{S}}&= 0, \end{aligned}$$

(45)

From the above equations, subgrid-scale kinetic energy can be written as

$$\begin{aligned} k_{\mathrm{SGS}}&= \frac{C_k}{C_{\epsilon }}\varDelta ^2(2{\bar{S}}:{\bar{S}}), \end{aligned}$$

(46)

$$\begin{aligned} k_{\mathrm{SGS}}&= \frac{C_k}{C_{\epsilon }}\varDelta ^2|{\bar{S}}|^2, \end{aligned}$$

(47)

where $|{\bar{S}}| = \sqrt{2 {\bar{S}}_{ij} {\bar{S}}_{ij}}$. Furthermore, substituting Eq. 40 in the above equation, we get

$$\begin{aligned} \nu _\mathrm{e} = C_k \varDelta ^2 \sqrt{\frac{C_k}{C_\epsilon }}|{\bar{S}}|. \end{aligned}$$

(48)

We can define a new constant coefficient as

$$\begin{aligned} C_\mathrm{s}^2 = C_k \sqrt{\frac{C_k}{C_\epsilon }}. \end{aligned}$$

(49)

where $C_\mathrm{s}=0.1678$ is called the Smagorinsky coefficient. Finally, we get following expression for SGS eddy viscosity

$$\begin{aligned} \nu _\mathrm{e} = C_\mathrm{s}^2 \varDelta ^2 |{\bar{S}}|, \end{aligned}$$

(50)

and the Smagorinsky model, given by Eq. 36, reads as

$$\begin{aligned} \tau _{ij}^{d} = -2C_\mathrm{s}^2 \varDelta ^2 |{\bar{S}}|{\bar{S}}_{ij}. \end{aligned}$$

(51)

Appendix B: Hyperparameters optimization

In appendix, we outline the procedure we followed for selection of hyperparameters for ANN with point-to-point mapping and neighboring stencil mapping. For ANN, there are many hyperparameters such as number of neurons, number of hidden layers, loss function, optimization algorithm, activation function, and batch size, etc. If we use regularization, dropout, or weight decay to avoid overfitting, the design space of hyperparameters increases further.

We focus on three main hyperparameters of ANN: number of neurons, number of hidden layers, and learning rate of optimization algorithm. The training data are scaled between $[-1,1]$ using the minimum and maximum value in the training dataset. We use ReLU activation function given by $\zeta (\chi ) = \text {max}(0,\chi )$, where $\zeta$ is the activation function, and $\chi$ is the input to the node. We use Adam optimization algorithm [71], and the batch size is kept constant at 256. Adam optimization algorithm has three hyperparameters: learning rate $\alpha$, first moment decay rate $\beta _1$, and second moment decay rate $\beta _2$. We test our ANN for two learning rates $\alpha =0.001$ and 0.0001. The other two hyperparameters in Adam optimization algorithm are $\beta _1=0.9$ and $\beta _2=0.999$. We employ mean-squared error as the loss functions, since it is a regression problem. We test both ANN with point-to-point mapping and neighboring stencil mapping for four different number of hidden layers $L=2,3,5,7$. The ANN with point-to-point mapping is tested for four different number of neurons $N=20,30,40,50$, and the local stencil mapping is tested for $N=40,60,80,100$. The number of neurons is higher in case of local stencil mapping because there are more features compared to point-to-point mapping.

The optimal ANN architecture is selected using multi-dimensional gridsearch algorithm coupled with k-fold cross-validation. Cross-validation is a procedure used to determine the performance of the neural network on unseen data. The procedure consists of dividing the training data into k groups, training the ANN by excluding each group and evaluating the model’s performance on that group. Therefore, if we use fivefold cross-validation, then the model is trained five times and the performance index is computed for five groups. Once the performance for each group is available, the mean of the performance index is utilized to select optimal hyperparameters. We use 500 epochs for determining the optimal hyperparameters. A good learning is achieved when both training loss and validation loss reduce till the learning rate is minimal. We apply coefficient of determination $r^2$ as the performance index to decide optimal hyperparameters. The calculation of coefficient of determination is done using the following formula

$$\begin{aligned} r^2 = 1 - \frac{\sum _{i}(y_i-{\tilde{y}}_i)^2}{\sum _{i}(y_i-{\bar{y}})^2}, \end{aligned}$$

(52)

where $y_i$ is the true label, ${\tilde{y}}$ is the predicated label, and ${\bar{y}}$ is the mean of true labels.

Figure 20 displays the performance index for ANN with point-to-point mapping and ${\mathbb {M}}3$ model for all hyperparameters tested using gridsearch algorithm. It can be observed that the performance of the network does not change significantly with hyperparameters and the difference in performance is very small. The optimal hyperparameters obtained for point-to-point mapping ANN are $L=2$, $N=40$, and $\alpha =0.0001$. We use the same hyperparameters for other two models ${\mathbb {M}}1$ and ${\mathbb {M}}2$ for point-to-point mapping ANN. We see the similar behavior in case of neighboring stencil mapping ANN and model ${\mathbb {M}}3$ as shown in Fig. 21. The optimal hyperparameters for neighboring stencil mapping ANN are $L=2$, $N=40$, and $\alpha =0.001$.

Fig. 20

Hyperparameters search using the gridsearch algorithm combined with fivefold cross-validation for the neural network using point-to-point mapping with ${\mathbb {M}}3$

Fig. 21

Hyperparameters search using the gridsearch algorithm combined with fivefold cross-validation for the neural network using neighboring stencil mapping with ${\mathbb {M}}3$

As discussed in Sect. 4.1, we get poor prediction between true and predicted stresses for point-to-point mapping with model ${\mathbb {M}}1$. Figure 22 shows the PDF of true and predicted stresses computed with different activation functions. It can be observed that the predicted stresses are almost the same for all activation functions. Therefore, we can conclude that we need additional input features such as velocity gradients to improve the prediction with point-to-point mapping.

Fig. 22

Probability density function for SGS stress distribution with point-to-point mapping. The ANN is trained using ${\mathbb {M}}1{:}\,\,\{{{\bar{u}},{\bar{v}}}\} \rightarrow \{{\tilde{\tau }}_{11},{\tilde{\tau }}_{12},{\tilde{\tau }}_{22}\}$ with different activation functions. The training set consists of 70 time snapshots from time $t=0.0$ to $t=3.5$, and the model is tested for 400th snapshot at $t=4.0$

The CNN architecture has similar hyperparameters as the ANN. Additionally, we need to select the kernel shape and strides for CNNs. Stride is the amount by which the kernel should shift as it convolves around the volume. We use the stride = 1 in both x and y directions. We use $3 \times 3$-shaped kernel in our CNN architecture. We check the performance of CNN architecture for different number of hidden layers $L=2,4,6,8$, different number of filters $N=8,16,24,32$, and two learning rates. Figure 23 displays the performance index of CNN for different hyperparameters. The performance of CNN is more sensitive to the learning rate, and we observe stable performance for the learning rate $\alpha =0.001$. The performance is almost similar for $L=6,8,10$ with different number of kernels. We can select $L=6$ and $N=16$, which has performance index of 0.76. Additionally, we test the CNN architecture with $L=6$ and [16, 8, 8, 8, 8, 16] distribution for the number of kernels along hidden layers and we observed the performance index of 0.75 at less computational cost. Therefore, we apply $L=6$, $N=[16,8,8,8,8,16]$, and $\alpha =0.001$ as our hyperparameters for the CNN architecture.

Fig. 23

Hyperparameters search using the gridsearch algorithm combined with fivefold cross-validation for CNN mapping with model ${\mathbb {M}}3$

Appendix C: CPU time measurements

In this study, the pseudo-spectral solver used for DNS is written in Python programming language. The code for coarsening of variables from fine to coarse grid, dynamic Smagorinsky model code is all written in Python. We use vectorization to get faster computational performance. The machine learning library Keras is also available in Python and is used for developing all data-driven closure models. Therefore, the CPU time reported in our analysis is for codes, which are all developed on the same platform. We would like to highlight that when the trained model is deployed, it makes the function for first time and hence it takes slightly more time. Once the function is created, the CPU time for deployment is less. Therefore, in all our tables, we report the CPU time for running the predict function second time since initializing CUDA kernels might yield a startup overhead as shown in Listing 1, where t1 here has some idle time due to initializing kernels. In our study, we report t2, and we further verified that t3 − t2 = t2, which illustrate that the reported CPU times are consistent.

Appendix D: ANN and CNN architectures

We use open-source Keras library to build our neural networks. It uses TensorFlow at the backend. Keras is widely used for fast prototyping, advanced research, and production due to its simplicity and faster learning rate. Keras library provides different options for optimizers, neural network architectures, activation functions, regularization, dropout, etc. Any simple neural network architecture can be coded with few lines of code. The sample code for ANN and CNN used in this work is listed in Listings 2 and 3.