2 Artificial neural networks

2.1 Architecture

The basic structure of the simple feed-forward ANN consists of layers possessing a predefined number of unit cells called neurons. Each of these layers has an associated transfer function and each unit cell has an associated bias. Any input to the neuron has a bias added to it followed by activation through the transfer function. To describe this process using matrix operations, we have for a single neuron in the $l$ th layer receiving a vector of inputs $\unicode[STIX]{x1D64E}^{l}$ from the $(l-1)$ th layer given by (Demuth et al. Reference Demuth, Beale, De Jess and Hagan2014)

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{l}=\unicode[STIX]{x1D652}^{l}\unicode[STIX]{x1D653}^{l-1},\end{eqnarray}$$

where $\unicode[STIX]{x1D652}^{l}$ stands for a matrix of weights linking the $l-1$ and $l$ layers with $\unicode[STIX]{x1D653}^{l-1}$ being the output of the $(l-1)$ th layer. The output of the $l$ th layer is now given by

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D653}^{l}=f(\unicode[STIX]{x1D64E}^{l}+\unicode[STIX]{x1D63D}^{l}),\end{eqnarray}$$

where $\unicode[STIX]{x1D63D}^{l}$ is the vector of biasing parameters for the $l$ th layer. Every node (or unit cell) has an associated transfer function which acts on its input and bias to produce an output which is fed forward through the network. The nodes which take the raw input value of our training data set (i.e., the nodes of the first layer in the network) perform no computation as they do not have any biasing or activation through a transfer function. The output layers generally have a linear activation function with a bias which implies a simple summation of inputs incident to a unit cell with its associated bias. In this investigation, we have used one hidden layer of neurons between the set of inputs and targets with a Tan-Sigmoid activation function which can be expressed as

(2.3) $$\begin{eqnarray}f(a)=\frac{1-e^{-2a}}{1+e^{-2a}},\end{eqnarray}$$

where the transfer function $f$ calculates the neuron’s output given its net input. In theory, any differentiable function can qualify as an activation function (Zhang, Patuwo & Hu Reference Zhang, Patuwo and Hu1998); however, only a small number of functions which are bounded, monotonically increasing and differentiable are used for this purpose.

2.2 Extreme learning machine

In this section we detail the extreme learning machine (ELM) approach to generalized single-layer feed-forward ANN training. This methodology was proposed in Huang, Zhu & Siew (Reference Huang, Zhu and Siew2004) for extremely fast training of a single-layer feed-forward ANN based on the principles of the least-squares approximation. We would like to note that the ELM training approach has previously been utilized for image deblurring but using the principles of classification in a classical deep convolutional network as against the principles of regression utilized here (Wang et al. Reference Wang, Huang, Luo, Wang and Luo2011). For the ease of description, let us define a few matrices for the single-layer feed-forward network. This is a generalization of the architecture introduced in § 2.1. Our input matrix is

(2.4) $$\begin{eqnarray}\unicode[STIX]{x1D653}^{0}=\left[\begin{array}{@{}cccc@{}}\boldsymbol{x}_{1}^{0} & \boldsymbol{x}_{2}^{0} & \ldots & \boldsymbol{x}_{n_{s}}^{0}\end{array}\right],\end{eqnarray}$$

where $\boldsymbol{x}_{i}^{0}$ is the $i$ th sample (out of a total of $n_{s}$ samples) of a multidimensional input column vector. Our weights connecting the inputs to the middle (hidden) layer and the biases associated with the hidden layer are given by

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D652}^{1}=\left[\begin{array}{@{}c@{}}\boldsymbol{w}_{1}^{1}\\ \boldsymbol{ w}_{2}^{1}\\ \vdots \\ \boldsymbol{ w}_{Q}^{1}\end{array}\right],\quad \unicode[STIX]{x1D63D}^{1}=\left[\begin{array}{@{}c@{}}\boldsymbol{b}_{1}^{1}\\ \boldsymbol{ b}_{2}^{1}\\ \vdots \\ \boldsymbol{ b}_{Q}^{1}\end{array}\right],\end{eqnarray}$$

where $\boldsymbol{w}_{i}^{1}$ and $\boldsymbol{b}_{i}^{1}$ represent row vectors corresponding to the weights and biases related to the $i$ th neuron (out of a total of $Q$ neurons) in the hidden layer. These are initialized to be small non-zero random numbers to enforce generalization. It is important to remark here that the extreme learning machine methodology prescribes biases only for hidden layer neurons (i.e., the output layer biases $\unicode[STIX]{x1D63D}^{2}=0$ ). The output of the hidden layer neurons becomes

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D643}^{\intercal }=f(\unicode[STIX]{x1D652}^{1}\unicode[STIX]{x1D653}^{0}+\unicode[STIX]{x1D63D}^{1}),\end{eqnarray}$$

where $f(\unicode[STIX]{x1D653})$ implies a Tan-Sigmoid activation procedure on each element of a matrix $\unicode[STIX]{x1D653}$ and the superscript $\mathsf{T}$ implies a matrix transpose. The weight matrix of the second layer may be given as

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D652}^{2}=\left[\begin{array}{@{}cccc@{}}\boldsymbol{w}_{1}^{2} & \boldsymbol{w}_{2}^{2} & \ldots & \boldsymbol{w}_{Q}^{2}\end{array}\right],\end{eqnarray}$$

where $\boldsymbol{w}_{i}^{2}$ now refers to the column vectors related to the hidden layer neurons. Our outputs of the ELM may thus be represented as

(2.8) $$\begin{eqnarray}\unicode[STIX]{x1D64E}^{2}=\unicode[STIX]{x1D652}^{2}\unicode[STIX]{x1D643}^{\intercal },\end{eqnarray}$$

as there is a linear activation in the outer layer with no bias. This output must be trained against a set of targets corresponding to each input vector given by

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D64F}=\left[\begin{array}{@{}cccc@{}}\boldsymbol{t}_{1} & \boldsymbol{t}_{2} & \ldots & \boldsymbol{t}_{n_{s}}\end{array}\right],\end{eqnarray}$$

where $\boldsymbol{t}_{i}$ is a column vector of the targets corresponding to the $i$ th sample. The ELM training mechanism is given as follows. In order to calculate the matrix $\unicode[STIX]{x1D652}^{2}$ , we must recognize that its optimal solution should satisfy

(2.10) $$\begin{eqnarray}\unicode[STIX]{x1D652}_{op}^{2}\unicode[STIX]{x1D643}^{\intercal }=\unicode[STIX]{x1D64F},\end{eqnarray}$$

or by taking a transpose of both sides

(2.11) $$\begin{eqnarray}\unicode[STIX]{x1D643}\unicode[STIX]{x1D652}_{op}^{2^{\intercal }}=\unicode[STIX]{x1D64F}^{\intercal },\end{eqnarray}$$

which leads us to the following expression for the optimal weights

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D652}_{op}^{2^{\intercal }}=\unicode[STIX]{x1D643}^{\dagger }\unicode[STIX]{x1D64F}^{\intercal }.\end{eqnarray}$$

The weights $\unicode[STIX]{x1D652}^{1}$ and biases $\unicode[STIX]{x1D63D}^{1}$ are generated initially using small random numbers. The utilization of small random numbers promotes generalization and prevents overfitting (Demuth et al. Reference Demuth, Beale, De Jess and Hagan2014). The matrix given by $\unicode[STIX]{x1D643}^{\dagger }$ is calculated using a generalized Moore–Penrose pseudoinverse (Albert Reference Albert1972; Serre Reference Serre2002). Once the optimal weights of the second layer are obtained, our network is trained for deployment. Due to the random number values chosen for the weights in the first layer, our network is well suited to highly effective generalization of the training process due to a smaller degree of freedom of the overall ANN. The proposed direct training procedure substitutes the gradient-based iterative optimization algorithms such as the Bayesian regularization method (MacKay Reference Mackay1992; Foresee & Hagan Reference Foresee and Hagan1997) in favour of the calculation of a pseudoinverse for the least-squares solution of optimal weights. This gives us the ability to obtain generalized predictive networks through training times which are several orders of magnitude lower in comparison with standard iterative techniques (Huang, Zhu & Siew Reference Huang, Zhu and Siew2006). Algorithms 1 and 2 describe the training and deployment process for the single-layer ANN trained using ELM. Effectively, Step 4 in Algorithm 1 represents the computational expense related to our supervised training.

3 Blind deconvolution

The primary motivation for the development of the aforementioned training and deployment procedures for our ANN architecture is described in this section. We are motivated by the approximate deconvolution framework (Stolz & Adams Reference Stolz and Adams1999), which uses an iterative resubstitution procedure known as the Van Cittert iterations for reconstructing the contribution of subfilter-scale content. This procedure may be represented tensorially as:

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ij}=\bar{u}_{i}\bar{u}_{j}-\overline{u_{i}u_{j}},\end{eqnarray}$$

where we obtain an approximation of these subfilter stresses as (Germano Reference Germano2015)

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D70F}_{ij}=\bar{u}_{i}\bar{u}_{j}-\overline{u_{i}^{\ast }u_{j}^{\ast }},\end{eqnarray}$$

in which the asterisk superscript indicates an approximately deconvolved variable. An important caveat of this mechanism is a user-defined low-pass spatial filtering kernel. While the general approach is to implement various versions of Gaussian blur kernel, it is quite possible that complex flow configurations may actually exhibit natural filter shapes which are much more contorted. Fortunately, the vast wealth of data that may be collected by modern experimental techniques represents an opportunity for us to bypass this user-defined heuristic stage. In essence, blind deconvolution refers to the estimation of the underlying blur kernel without any knowledge or closed form estimate of its true nature. We must emphasize here that while there exist several algorithms that are truly blind in nature (for instance in Dabov et al. (Reference Dabov, Foi, Katkovnik and Egiazarian2007)), this investigation proposes a framework which is purely data-driven. To summarize, our data may be considered to be an aid in implicitly estimating the shape of the blur kernel through the form of a trained neural network.

Our supervised learning framework requires training data which encapsulates the relationships between the coarsened true data and its perturbed versions. Note that, as is common in data-driven strategies, we normalize all our data to a range between $-1$ and 1. This choice of the normalized data range is due to the fact that the Tan-Sigmoid activation function provides outputs between these limits as well. The perturbation may be introduced through a low-pass spatial filtering procedure:

(3.3) $$\begin{eqnarray}u_{i}^{\prime }=G(\boldsymbol{x};\unicode[STIX]{x1D70E})u_{i},\end{eqnarray}$$

which utilizes the following Gaussian kernel with standard deviation $\unicode[STIX]{x1D70E}$ (held at a default value of 1.0):

(3.4) $$\begin{eqnarray}G(\boldsymbol{x};\unicode[STIX]{x1D70E})=\frac{1}{(\sqrt{2\unicode[STIX]{x03C0}}\unicode[STIX]{x1D70E})^{d}}\exp \left(-\frac{\left|\boldsymbol{x}\right|^{2}}{2\unicode[STIX]{x1D70E}^{2}}\right),\end{eqnarray}$$

where $d$ is the number of dimensions of the field being blurred. The use of a low-pass spatially filtered field tests the ability of the proposed ANN architecture for inertial range recovery and becomes a true test of the blind deconvolution ability of a neural network.

Another perturbation which is examined in this investigation is the addition of a normally distributed noise with an amplitude on the order of a tenth of the maximum value of the field variable which can be expressed as:

(3.5) $$\begin{eqnarray}u_{i}^{\prime }=u_{i}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D705},\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ is a normally distributed random number between $-1$ and 1 and the coefficient $\unicode[STIX]{x1D707}$ (with a default value of 0.2) corresponds to the magnitude of the noise added to the field. The purpose of training relationships between inputs with noisy perturbations and the true solution is to test the stability of this approach for high-wavenumber energy accumulation. This is particularly important for coarse-grained large eddy simulations which are susceptible to numerical overflow. Thus, our supervised learning framework utilizes $u_{i}^{\prime }$ as our training inputs and $u_{i}$ as our targets. Note that we train separate networks for each type of perturbation. The ELM training mechanism leads to tractable training times even for cases with several hundred neurons in the hidden layer. The results obtained in this study have employed 100 neurons for training and testing data obtained from different realizations of field values. A schematic of the blind deconvolution methodology integrated into the one-layer ELM-training-based ANN is shown for a two-dimensional test case in figure 1.

Figure 1.

A schematic of the proposed blind deconvolution strategy for a two-dimensional test case where we aim to recover $\unicode[STIX]{x1D714}^{\ast }$ from $\bar{\unicode[STIX]{x1D714}}$ . Note that each hidden neuron is associated with a unique bias.

4 Cross-validation

An important assessment of a data-driven modelling approach is a cross-validation procedure to estimate the performance of the model for situations which it has not been trained for. In other words, it is necessary to ascertain that the data-driven model is not localized to the training data and provides a similar performance for different solution fields exhibiting similar physics but different numerical arrangements or magnitudes. Therefore, it is vital to test the performance of the trained ANN on a set of data which it has not seen previously (during training). For the purpose of data generation, we utilize high-fidelity simulations for our previously mentioned test cases so that their flow physics may be accurately resolved. Our 3D numerical experiments are generated using $512^{3}$ degrees of freedom while our 2D case is generated using $2048^{2}$ grid points. A coarse-grained large eddy simulation is mimicked through the subsequent subsampling of these high-quality data sets. Our coarse-graining procedure leaves us with a subsampled field of $64^{3}$ and $256^{2}$ degrees of freedom in the 3D and 2D test cases, respectively. Therefore, our coarse-graining procedure involves the selection of every eighth grid point in the high-fidelity uniform grid data.

We may devise several shifted data sets (each possessing $64^{3}$ degrees of freedom in 3D and $256^{2}$ degrees of freedom in 2D) but which continue to represent the physics of the fine-grained dataset. Missing data points at boundaries may be reconstructed through the use of the periodic boundary conditions for the given test cases. A simplified shifting schematic is shown in figure 2, where a simple technique of generating two coarse data sets from fine data is demonstrated. In a nutshell, the coarse-graining and shifting procedures allow us to devise up to 63 and 511 completely different data sets in two and three dimensions, respectively. For the purpose of cross-validation we randomly choose any four of these multiple data sets for the generation of three sets of testing data and one set of training data.

Figure 2.

A schematic of the spatial shifting strategy for a simplified two-dimensional grid showing two different data sets. Several different data sets may be constructed using this technique which are numerically different physically similar.

Once our four different randomly generated data sets are identified, we add perturbations corresponding to the type of behaviour demanded of the proposed artificial neural network. To test the deconvolution ability of our network, our training data (i.e., one of the four data sets) is filtered with an appropriate Gaussian smoothing (for our inputs to the network) and unfiltered training data are utilized as outputs to the network. Three testing data sets are generated in a similar manner, with one of these data sets being filtered with the same filter radius and the others being filtered with a 10 % larger and 10 % smaller filter radius, respectively. The trained ANN is then utilized to recover deconvolved approximations to the true field for these three test data sets. This ensures that the trends of the trained network are not due to overfitting or ‘data-memory’ but through an implicit learning of the inverse filtering. A similar procedure is also utilized to cross-validate the regularization ability of the closure wherein the randomly chosen data sets are perturbed through Gaussian noise. A table describing the filter radii and magnitudes of noise for our training and testing data sets is shown in table 1. Note that these data sets are all generated from the same high-fidelity solution field and correspond to a perfectly a priori analysis. A concise summary of this cross-validation procedure is given in a flow chart in figure 3.

Figure 3.

A flow chart detailing the generation of cross-validation data sets using the shifting operation as well as different perturbations. The network is trained using Algorithm 1 and is then tested on test data 1, 2 and 3.

Table 1.

Cross-validation data sets for the proposed data-driven blind deconvolution closure.

In addition to the a priori cross-validation outlined above, we also assess our proposed architecture by utilizing training and testing data across different flow physics. This procedure is carried out to test the universal nature of the learning for the purpose of both deconvolution and regularization abilities. This may be assumed to be another challenging tier for cross-validation. We elaborate this step in further detail in § 6.

6 Universality

The true value of the implementation of a data-driven framework emerges from an intelligent combination of our knowledge of physics from first principles with the robust techniques inherent to data-driven modelling. One of the motivations of studying homogeneous isotropic turbulence test cases in this investigation was to link the universal nature of the physics (as evidenced in the $k^{-5/3}$ scaling observed in the averaged kinetic energy spectra) through our proposed framework across flows exhibiting distinctly different physics. For instance, a great utility of any data-driven modelling framework would be to utilize training data from lower Reynolds numbers or shorter (high-fidelity) simulations to devise data-driven closures for runs with higher Reynolds numbers or longer simulations which are computationally intractable.

We first examine the ability of the closure modelling strategy by using training data from a high-fidelity Taylor–Green vortex simulation at $Re=1600$ to recover true fields at $Re=5000$ as shown in figure 19. We observe that the framework is successfully able to replicate closure performance as shown in previous test cases. It tells us that the deconvolution (or regularization) ability has been learnt without violating (as it appears) the underlying physics of homogeneous isotropic turbulence. Our next assessment shown in figure 20 utilizes the closure modelling performance of the framework when training data is utilized at time $t=15$ for a Taylor–Green vortex test case to obtain true field reconstructions at $t=20$ (for the same simulation). It can once again be observed that the closure recovery is exceptional, with trends similar to those exhibited by previous test cases.

Next, we present cross-validation results for different flows. Here, the proposed blind deconvolution ANN closure is trained from a completely different simulation data (i.e., testing with different flow configurations). Figure 21 demonstrates the ability of the data-driven closure to predict subfilter terms for the stratified KHI turbulence test case where TGV simulation data at $Re=1600$ is utilized to train our ANN. In this case we highlight that the only link between training and testing data is the $k^{-5/3}$ cascade and the magnitude of perturbations to the coarse-grained field. It can be seen that a remarkable inertial range recovery and regularization performance is still observed for both the low-pass spatially filtered and noisy testing data. Analogously, in figure 22, we invert our training and testing data with stratified turbulence data being used to train the ELM, which is further used for deconvolution and regularization of an appropriately perturbed solution field given by a Taylor–Green vortex simulation at $Re=1600$ . Once again, previously observed trends in the output of the closure model are recovered. From the evidence provided above, we may tentatively conclude that for flows exhibiting similar cascade characteristics, subfilter recovery and high-wavenumber regularization are reliant on the filter radius and magnitude of high-frequency perturbation alone. This allows for exciting possibilities for 3D turbulence closures, for instance, to implicitly learn filter shapes from physical experiments and reproduce accurate trends through coarse-grained numerical simulations within a wide range of Reynolds numbers.

7 Concluding remarks

An artificial neural network architecture is proposed for the data-driven deconvolution and regularization of low-pass spatially filtered turbulence fields. Both 2D and 3D test cases are examined with training data sets obtained through the coarsened and perturbed versions of high-fidelity simulations for canonical homogeneous isotropic turbulence and stratified compressible turbulence problems. Two types of perturbations are tested, one in which the true data is filtered using a Gaussian kernel and the other which has a certain quantum of noise added to the field to represent high-wavenumber aliasing errors. Testing data sets are generated through a spatial shifting procedure as well as through the utilization of slightly different magnitudes of filter radius and noise. This ensures that deconvolution and regularization performance of the proposed architecture remains localized to physics and not numerical artefact. The proposed architecture is tested independently on Kraichnan, Kolmogorov and stratified compressible turbulence test cases and is able to provide an estimate of the deconvolved variable as examined through the energy spectra and probability density functions of the recovered fields. For noisy data, the innate regularization of the ELM training approach results in smoother predictions for the recovered variable. For our training mechanism, the ELM approach is chosen due to its exceptional speed of training in comparison with traditional gradient-based methods as well as its excellent generalization ability (which avoids the problem of overfitting). The ability of the proposed approach for subfilter-scale content recovery is benchmarked against several popular structural closure modelling strategies. It is observed that our the data-driven framework yields a similar deconvolution performance without the explicit specification of a filter kernel. In addition, we also test the universality of our data-driven closure by utilizing testing and training data sets from different simulations displaying the same classical Kolmogorov cascade. This leads to an exciting observation, as described in § 6, which tells us that the proposed framework may be utilized to leverage our understanding of the cascade of energy in three-dimensional turbulence to obtain reconstructions for flows across different physics. Indeed, the blind deconvolution procedure appears to be linked across a wide range of physics solely through the underlying filter radius that arises in LES flow computations due to coarse-graining and implicit (or explicit) numerical dissipation.

A natural follow-up to this investigation is to test our proposed approach in a fully a posteriori analysis. One of our primary goals in subsequent investigations is also to address the issue of sampling for training data. Sampling strategies must be devised to ensure that the proposed framework is exposed to data from many physical regimes (as opposed to indiscriminate selections in parameter regimes). This is an important distinction that has significant implications on the performance of any proposed data-driven modelling framework. Basically, a data-driven model is only as good as the data it has been trained on and can only reproduce physical behaviours similar to those it has seen in training. In addition to sampling strategies, it is also important to develop outlier identification systems for noisy data. Data preprocessing, an active area of research in the data science community, must therefore be integrated into our framework. From this point of view, the fast training times of the proposed architecture suggest the use of multiple networks trained to act in a ‘committee’ for aggregate subfilter predictions of flow datasets with variable spatial and temporal characteristics. This ensures robustness towards outliers, since aggregate predictions are likely to be more accurate than those from solely one trained network.

Since, apart from its training data, an ANN’s performance is heavily dependant on its architecture, it is necessary to estimate the performance of the proposed framework for several different configurations in terms of network architectures (number of inputs, number of outputs, number of neurons, activation functions etc). Fortunately, the speed of training aids us in generating regularized network weights for extremely large datasets with considerable ease. In conclusion, this work features a preliminary glimpse at an exciting avenue for the next generation in data-driven turbulence closures.