Turbulence Model Development based on a Novel Method Combining Gene Expression Programming with an Artificial Neural Network

GEP [12] is a powerful EA which can solve a series of problems, such as symbolic regression, sequence induction with and without constant creation, block stacking, etc. For the symbolic regression task, a population of mathematical expressions is evolved through natural selection and genetic operations over many generations. Unlike traditional genetic algorithms, the fundamental advantage of GEP lies in its unique design of the individual, i.e., the chromosome, which contains mathematical expression information. Each chromosome has a distinct genotype and phenotype. The purpose of the genotype is to store genetic information that can be evolved in each generation, while the phenotype corresponds to the explicit expression extracted from the genetic information.

The genotype of the chromosome is encoded as fixed-length strings, which consist of one or more genes and a set of linking symbols to link genes. A single gene can represent a mathematical expression with a string of symbols. For example, a polynomial of $x$ and $y$

can be encoded as a single gene with a genotype string

where the mathematical symbols $\times$, $+$, and $-$ represent multiplication, addition, and subtraction, respectively, the alphabetic characters $x$ and $y$ represent the variables, and $1$ and $2$ represent constant numbers.

The genotype string in GEP can be decoded into its phenotype, i.e. a mathematical expression, with an expression tree (ET) structure. As shown in Fig. 1, the genotype string given by Eqn. (2) can be decoded by taking the root symbol, filling its child nodes with the following $n$ symbols, and repeating until no more symbols can be added to this tree. Here, $n$ represents the arity of the symbol, e.g., multiplication $\times$ takes two arguments, so $n_{\times}=2$.

By combining multiple genes, the chromosome can represent long and complex expressions. It can obtain its phenotype by decoding its genes separately and connecting those ETs with link symbols. Consider a chromosome with two genes $<gene1>$, $<gene2>$ and a link symbol $\times$. The genotype of this chromosome can be written as

The phenotype can be expressed as shown in Fig. 2. An expression described by the ET can be calculated directly from the leaf nodes to the root node, which is easy to implement in the program.

Due to the design of the genotype, a chromosome can easily participate in genetic operations because all the genetic operations applied to a chromosome are equivalent to replacing part of the gene with another fixed-length string. When this chromosome needs to be evaluated, it can be recursively decoded to the phenotype and efficiently calculated as a regular mathematical expression.

Compared to other data-driven algorithms, for example, ANN [3, 23] and random forest [24], the advantage of GEP is its global searching strategy. GEP can generate complex mathematical expressions by applying simple genetic operations, and the variation in chromosomes allows exploring the entire search space. Thus, GEP is a global optimization algorithm that can search the optimization space and has the potential to address the problem of function discovery and symbolic regression. However, genetic operations do not consider the optimization direction from the training data, leading to convergence difficulties. Moreover, constants in the expression depend on the combination of preset symbols, making it difficult to obtain accurate values of the constant numbers [25, 26].

SRNN [5] is an ANN architecture recently designed for symbolic regression tasks. Fig. 3 exhibits the schematic of an example SRNN, with two hidden layers for visual simplicity. The basic structure of SRNN is usually a fully connected ANN, in which the $i^{th}$ layer of the network can be derived from the $(i-1)^{th}$ layer:

where $I_{i-1}$ is the input data from the $(i-1)^{th}$ layer, $W_{i}$ is the weight matrix of this connection, and $O_{i}$ is the output data of the $i^{th}$ layer. It is noted that in SRNN, the activation functions $f(*)$, which are typically nonlinear mappings like sigmoid or Relu [27] in traditional ANNs, are simple mathematical operators. Specifically, the activation function $f(*)$ in the layer $i$ may consist of separate functions for each node of layer $i-1$ (such as $sin$ or $ln$) and may include functions that can adopt more than one argument to generate one output (such as the multiplication function).

SRNN can fit complex combinations and compositions of the input variables with various primitive functions by stacking multiple layers and applying different activation functions. It can be trained end-to-end with backpropagation and other powerful deep learning techniques while still producing interpretable expressions. However, gradient-based methods suffer from a gradient-loss problem at the local optimum or saddle point during training. Moreover, the network structure, which can significantly affect the performance of the training, needs to be decided before training and is not flexible, thus limiting the application of this approach [28, 29].

In order to show the advantage and disadvantage of the GEP and SRNN methods, we test their performance for a simple polynomial case. Consider the following polynomial expression:

where $x_{1}$ and $x_{2}$ are the input variables, and $C_{1,\dots,5}$ are constants. By setting the constants as

Eqn. (5) can be expanded and simplified as follows:

which is used to generate a training dataset with $x_{1},x_{2}\sim N(0,1)$. Note that random noise with a magnitude of $10\%$ is then superposed to the outputs y.

For GEP, the chromosome of the expression in Eqn. (5) contains five genes, and the genotype can be written as

The initial population for the GEP training is randomly generated, except for that an individual with chromosome Eqn. (LABEL:eq:poly_gene) is created and inserted into the initial population, which ensures that the starting point of GEP and SRNN is consistent. The values of $C_{1,\dots,5}$ are initialized with a normal distribution $N(0,1)$. It is noteworthy that multiple genotypes can exist for the same expression.

Furthermore, two different SRNNs are built to reflect the influence of the ANN architecture on the performance, as presented in Fig. 4. SRNN1 is designed explicitly with the prior knowledge of Eqn. (5), and the mathematical operators are selected to ensure that SRNN1 can represent this polynomial. Thus, SRNN1 can easily represent Eqn. (5) by setting the special connection weights. By contrast, SRNN2 is much simpler and arbitrary, leading to a limited exploration space. The connection weights of SRNN1 and SRNN2 are also initialized with a normal distribution.

Fig. 5 compares the training processes of the GEP, SRNN1, and SRNN2 for Eqn. (5). The root-mean-squared error (RMSE) of SRNN1 decreases quickly during the training stage, and its training terminates after $1,000$ epochs, reaching a low error value. Owing to its simpler architecture, SRNN2 cannot identify the correct form of the target expression, leading to convergence to a higher RMSE, even after $30,000$ epochs of training. On the other hand, SRNN1 can extract constants from the corresponding connection weights

which are very close to the target values of Eqn. (6). This shows that using the gradient information of the training data, the SRNN can determine the optimization direction and adjust the constants accordingly. By setting a reasonable network structure, SRNN1 converges to the best result. However, with an improper SRNN structure, the SRNN2 shows much worse results.

As a comparison, the RMSE value given by the GEP training results in Fig. 5 is still quite high after $1,000$ generations, which costs much more computational resources compared to SRNN1. The expanded form of the GEP result is

which represents the same functional form as in Eqn. (7). However, the constant values of the GEP expressions are inaccurate because searching the global optimization space and generating precise constant numbers by traditional genetic operations are inefficient.

The results of the polynomial case reflect two aspects of the symbolic regression task: one is to uncover the underlying expression structure, which is an advantage of GEP; and the other is to generate a proper set of constant coefficients in the expression, which can be accomplished efficiently and accurately using the SRNN. At the same time, the presented results also suggest that the GEP cannot consider the optimization direction, leading to difficulty in obtaining accurate constants and thus slowing the overall convergence rate. The SRNN performance, on the other hand, is greatly affected by the preset network structure because the structure defines the search space and the most complex expression form that the SRNN can discover. However, there is usually a trade-off between model complexity and network performance, and the setting of the network structure usually requires prior knowledge of the application fields. Therefore, we intend to develop a novel method for general symbolic regression tasks in the present study by combining the advantages of GEP and SRNN, with the capability of efficiently discovering the expression structure and regressing the model coefficients at the same time.

The critical challenge for combining GEP with SRNN is to transform expressions between GEP and SRNN. An expression in GEP is encoded as the genotype of a chromosome, which is essentially an ordered string. However, in a SRNN, a structured network with weighted connections contains all the information of the expressions. To optimize the expression with the SRNN, we need to transform the genotype of chromosomes into ANN architectures. In other words, the SRNN should be initialized using GEP chromosomes. At the same time, to continue the GEP training iterations after the SRNN optimization step, we need to be able to transform the SRNN back into the GEP genotype without breaking the original structure of the chromosome. Therefore, the transformation of expressions between GEP and SRNN requires a special design for GEP individuals and the transformation principles.

Fig. 6 shows the design of GEP and SRNN architectures for transformation. Essentially, the transformation of the GEP into the SRNN is a process similar to decoding a chromosome into its ET, which requires a recursive analysis of each symbol. As an example, the expression

has the following genotype:

where $<*>$ represents a gene, $C_{i}$ are constant numbers, and $x_{i}$ denote the input variables. As discussed in Section 2.1, the genotype in Eqn. (12) can be decoded into the phenotype in Eqn. (11) via the ET shown in Fig. 6(a), which is then directly linked to the SRNN architecture with input neurons $x_{1}$, $x_{2}$ and $1$ shown in Fig. 6(b).

Specifically, the variable symbols $x_{i}$ in chromosomes are transformed into the input neurons in SRNN, and the operator symbols $\{+,-,\times\}$ are converted as hidden neurons and connected according to the sequence of chromosomes. Moreover, the constant symbols $C_{i}$ in GEP chromosomes are transformed into the connections between the input neuron $1$ and the hidden layers, indicated by the blue solid lines in Fig. 6 (b). In addition to the existing constant symbols in the GEP chromosome, each connection in SRNN, which is indicated by the red solid lines in Fig. 6 (b), contains a trainable weight. Thus, we can determine a special coefficient for each symbol, shown as red nodes in Fig. 6 (a). These coefficients are initialized with a value $1$ and can be optimized by SRNN.

Based on the principle listed above, the expressions can be converted from GEP to SRNN as shown by the corresponding ET and SRNN structures in Fig. 6. Thereafter, the resulting SRNN can be trained end-to-end using the gradient information obtained from the training data to optimize the trainable weights. Benefitting from the transformation of the GEP into the SRNN, the constant values in the mathematical expression of the chromosome can be promptly optimized to their local optimum. It should be noted that as traditional GEP has difficulty training constant coefficients. The number of constants in GEP models is usually limited by the preset constant symbols and the maximum length of chromosomes, and even inaccurate unless with exceptionally long training, as discussed in Section 2.3. Therefore, the introduction of coefficients for every symbol is usually not considered in traditional GEP. The extension of the symbol coefficients brings a stronger exploration ability for GEPNN, which is expected to significantly enhance the expressivity of chromosomes, especially the ability to train accurate model coefficients.

It is also noted that the structure of the hidden layers of the SRNN is identical to the green-boxed part of the ET in Fig. 6(a). Therefore, after the SRNN optimization, the transformation of the SRNN back to the GEP ET and hence the chromosome is straightforward. Subsequently, the expression optimized by SRNN can be included in the evolution of the next GEP generation, and new created expression structures generated by GEP can be further optimized based on the GEP-SRNN transformation in this iterative process. In this way, by combining GEP with SRNN, we aim to efficiently and accurately train the expression structure and the coefficients at the same time.

Based on the transformation process of expressions between GEP and SRNN discussed above, the framework of GEPNN is introduced in Fig. 7, in which the expression structures are determined using genetic algorithms, and the constants are efficiently optimized with SRNN. Initially, a population of individuals is randomly created. Each individual’s fitness is then evaluated according to the user-defined loss function, i.e. the training objective. If the termination criterion (e.g., achieving a preset generation number or obtaining a loss lower than a preset value) is not fulfilled, individuals compete with each other depending on their fitness, and the winner is selected for the mating pool individuals. Next, the individuals are updated by applying genetic operations to the mating pool. For the traditional GEP method, the updated population is evaluated, and the training process continues over the next generations. However, for the GEPNN, with the integration of the SRNN, the individuals are selected, transformed and optimized using ANNs. After SRNN optimization, these ANNs are recovered to chromosomes and updated the corresponding individuals to participate in the subsequent generations.

Traditional genetic operations in GEP include mutation and crossover. All individuals in the mating pool can perform one or more operations to manipulate and change their chromosomes. We remark that the SRNN optimization introduced in the GEPNN can be considered as an additional genetic operation specifically designed to enable efficient optimization of the constants in GEP chromosomes. Therefore, in the GEP iterations, all individuals in the mating pool have the chance to perform the genetic operation and SRNN optimization. In addition, $m$ candidates with the highest fitness values, i.e., those with the best expressions, are automatically subjected to SRNN optimization. Thus $m$ is a hyperparameter that depends on the population of the GEPNN and computational resources.

It is noteworthy that SRNN optimization is not applied in every generation of the GEPNN, as it incurs additional computational costs and is usually unnecessary. The SRNN optimization is applied for every $n$ generations, and $n$ is another hyperparameter that can be modified depending on the training situation, for example, training data size, and complexity of the training target. In addition, the SRNN optimization is applied when a new expression structure with better fitness values appears after other genetic operations.

In principle, GEPNN compensates for the disadvantages of GEP and SRNN. Fig. 8 shows a schematic of the optimization processes of the GEP, SRNN, and GEPNN. GEP uses genetic operations for mathematical expressions to create new forms of equations for global searching. However, there is no direction for optimization in the subsequent generation. Unlike GEP, SRNN uses a local gradient to instruct the equations to descend to the local optimum. However, it may fall into a local optimum or saddle point without gradients, and the network structure is critical to its performance as shown in Section 2.3. By integrating GEP with the SRNN, GEPNN combines the advantages of the two methods to obtain diverse individuals that can be optimized with gradients, finally finding the global optimum efficiently and accurately.

In addition to the GEPNN framework, various strategies have been introduced to enhance the sparsity of the trained models. We remark that enhancing model sparsity is usually desired for ML tasks, as it has the potential to make the trained model more interpretable and also avoid overfitting  [20]. For the GEPNN framework, ANNs are incorporated in the training iteration, so it is natural to apply existing strategies for the sparsity of ANNs here. The details of the sparsity strategies introduced to the present GEPNN method are thus discussed below.

First of all, for the training of SRNN, a regularization term is introduced to the loss function, and the intention is to set as many connection weights to $0$ as possible so that those connections can be removed from the final trained model expression. Adding a regularization term to the loss function is a straightforward and popular method to enforce the sparsity of neural networks [30, 31]. The commonly used $L_{q}$ regularization acts as:

where $w_{i}$ denotes the connection weights in the neural network. In the present study, we use a smoothed $L_{0.5}$ regularization [5] to avoid the singularity of the gradient as the weights go to $0$

Here, $a$ is the transition point between the standard $L_{0.5}$ regularization and the smoothed function, which is set to $a=0.01$ in our experiments.

Second, we use a multiphase training strategy following [5] and [32] for the SRNN optimization, which divides the training into two stages. In the first stage, a standard learning rate is set to optimize all the connection weights in the SRNN. After this stage, weights below a certain threshold $\beta$ are set and fixed as $0$, which filters the small-magnitude weights out and enforces the sparsity of SRNN. Then, in the second stage, the training of the filtered SRNN continues with a reduced learning rate to fine-tune the remaining weights. This multiphase training strategy ensures the sparsity of SRNN and allows the sparse network to be fully trained.

Finally, it is important to make sure that the sparsity enforced during the SRNN optimization can be kept after the expressions are transformed into GEP chromosomes. Therefore, we introduce non-coding symbols to prune the chromosome. The non-coding symbol is defined here as a symbol without expressive function. Similar to the GEP algorithm itself, the idea of non-coding symbols is also from the biology field, in which the non-coding DNA represents the sequences that do not encode proteins. A non-coding symbol in GEPNN will be created only during the expression transformation from SRNN to GEP. In this process, if a connection weight in SRNN is below a threshold $\beta$, its corresponding symbol in the chromosome will be replaced with a non-coding symbol, which means the coefficient of this symbol is small enough to be ignored. In particular, to maintain the structure of the chromosomes, the non-coding symbols are distinguished into two categories, i.e. the non-coding operators $N_{o}$ which have more than one arity and will take the place of the operator symbols, and the non-coding variables $N_{v}$ which will only appear on the leaf nodes. Fig .9 shows an example ET and the corresponding SRNN structure of a pruned expression. While calculating its fitness, the non-coding symbols return value $0$, which act as pruning. By introducing non-coding symbols, the sparsity of GEPNN is further strengthened. Note that the non-coding symbols can be replaced with valid symbols via genetic operations in GEP, which ensures that the pruned expression has the potential to evolve in the following generations.

With the regularization and non-coding symbols, GEPNN can obtain simpler and more interpretable expressions without affecting the training of SRNN and the evolution of GEP.

The GEP framework employed in this study was developed by Weatheritt and Sandberg [13], and its implementation was programmed using Python [33]. The SRNN used in this study is implemented via the Tensorflow, which is a powerful framework for deep learning. In the following paragraphs, we will illustrate the details and hyperparameters of GEPNN.

The loss function of the SRNN is defined using the mean-squared error between the prediction $\hat{y}$ and label $y$ using the mini-batch gradient descent algorithm, which can be written as

where $N$ denotes the batch size. The Adam optimizer [34], which considers the historical effects of gradients, is used to update the weights in the ANN.

For gradient-based approaches, the learning rate is an important hyperparameter. In the early training process of GEPNN, the learning rate of SRNN is set to a relatively large value to accelerate the optimization, while for later stages the learning rate of SRNN needs to be fine-tuned. Thus, the learning rate changes during the GEP training as $L_{r}^{(g)}={L_{r}}/{g}$, where $g$ refers to the $g^{th}$ generation of GEP iterations and $L_{r}$ is the global learning rate as a hyperparameter. Furthermore, for each GEP generation with fixed $g$, the learning rate of each epoch in the SRNN optimization exponentially decays as follows:

where $L_{r}^{*}$ and $g_{e}$ are the learning rate and global epoch number in the current training step, respectively. Moreover, $L_{r}^{(g)}$ is the initial learning rate in this generation, and $d_{r}$ and $d_{e}$ are hyperparameters representing the decaying rate and decaying epoch, respectively. It is noted that with the decaying learning rate, the training speed at the beginning is relatively fast, while numerical stability can be ensured approaching the final result.

The number of training epochs of each SRNN optimization step are adjusted dynamically according to convergence. After a fixed number of training epochs, we use the loss value to determine whether the training is converged. The training process of SRNN will be terminated if (1) the current loss is less than a preset threshold, (2) the loss value no longer decreases, or (3) the training epochs achieve the preset maximum epochs.

Moreover, the threshold $\beta$ used for masking the SRNN and pruning the expression in Section 3.3 is dynamically adjusted according to the magnitude of weights, defined as follows:

where $b$ is a hyperparameter denoting the threshold scale.

The settings of all hyperparameters are listed in Table 1. We set a large value for total generations to ensure that the GEP and GEPNN converges during the training stage. Owing to the different training expenses of genetic algorithms and ANNs, we use the CPU time consumed during training to compare the performance of each method.

In this subsection, we test the GEPNN method by training three physics laws based on synthetic data. The first case, the ideal gas law, is basic and consists of two input variables. The second case (i.e., the gravitation law) is more complex, involving three input variables and requiring an additional data-normalization operation. The third case, the logarithmic law in a turbulent channel flow includes various degrees of noise, testing the robustness of the GEPNN. Table 2 summarizes the setup and parameters of these three cases.

In this section, we apply the GEPNN to subgrid scale (SGS) stress modelling in large-eddy simulation (LES). LES is an important method for the numerical simulation of turbulence and has been applied widely in various applications in, for example, automobile design and aerospace engineering. By filtering the governing equations, the LES solves the large-scale flow field and models the effects of SGS structures on the large-scale flow field. These effects are considered by modelling the stress term $\tau_{ij}$ in the filtered governing equations, and this type of model is called SGS closure. Accurate modelling of the SGS stress is essential for turbulence simulation to predict the precise effects of subgrid-scale structures. With the development of LES, different models have been developed, including the Smagorinsky model [38], the dynamic Smagorinsky model (DSM) [39], the dynamic mixed model (DMM) [40], and the gradient model [41], etc. Although these models have been tested under different flow conditions and shown some success in predicting mean flow behavior, the physical quantities predicted by these models sometimes have low correlations with those of the DNS results. Using data-driven approaches to model the SGS stress aims to predict physical quantities more consistent with the DNS data.

The governing equations for 3-D incompressible turbulence [37, 42] are as follows:

where $u_{i}$ is the velocity component in the $i^{th}$ direction, $p$ is the pressure, $\nu$ is the kinematic viscosity, and $\mathcal{F}_{i}$ is a large-scale forcing.

In LES, the governing equations are filtered to separate large-scale motion and SGS stress. By applying the filter operation to Eqn. (27), the filtered governing equations [43] become

where the overbar variable $\widetilde{*}$ is referring to filtered quantities, and the SGS stress tensor $\tau_{ij}$ can be expressed as

This term is not closed because $\widetilde{u_{i}u_{j}}$ cannot be expressed with quantities from the filtered equations. Therefore, the SGS stress tensor $\tau_{ij}$ has to be modelled using known quantities from the flow field. According to the Cayley–Hamilton theory following Pope [44], the anisotropic part of the SGS stress tensor $\tau_{ij}^{A}$ can be expressed as

where $T^{(n)}$ is a tensor basis function, ${\widetilde{S}}_{ij}$ is the filtered strain-rate tensor, ${\widetilde{\Omega}}_{ij}$ is the rotation-rate tensor, and $I_{i}$ are independent invariants. The coefficients $g^{(n)}$ are functions of the invariants, which must be predicted using data-driven approaches.

The filtered strain rate tensor ${\widetilde{S}}_{ij}$ and the rotation rate tensor ${\widetilde{\Omega}}_{ij}$ are defined as follows:

Using the inverse timescale $|\widetilde{S}|$ to nondimensionalise ${\widetilde{S}}_{ij}$ and ${\widetilde{\Omega}}_{ij}$, the dimensionless quantities $s_{ij}$ and $\omega_{ij}$ are used to build the SGS model

The first four normalized tensor basis functions $T^{(1\dots 4)}$ and the first four normalized invariants $I_{1\dots 4}$ are selected:

where $I$ denotes an identity matrix. For further details on the training setup, refer to  [14].

To compare the differences between GEP, SRNN, and GEPNN, two SRNNs with different architectures, shown in Fig. 13(a) and Fig. 13(b), are designed and used for SGS stress modelling. Different from the SRNNs in Sec. 2.3, four separated SRNN blocks are used to predict the four basic coefficients $g^{(1\dots 4)}$. The SRNN blocks are connected with their tensor basis functions to generate the SGS stress and the gradient can be backpropagated to each block to ensure that all the blocks can be trained simultaneously.

The training dataset is extracted from the homogeneous isotropic turbulence [45, 14]. The grid size of the DNS is $1,024^{3}$, and the Taylor Reynolds number is $Re_{\lambda}\approx 260$ [46]. The filtered velocity, corresponding SGS stress $\tau_{ij}$, and other filtered variables are obtained using a top-hat filter with a width of $\mathrm{\Delta}=16\mathrm{\Delta}x$. Using all the data points in the filtered DNS (fDNS) data for training is inefficient and redundant, thus we downsampled the fDNS data and selected about $2\times 10^{6}$ pairs of data points for the training process.

During the training stage, the correlation coefficient $C(\tau)$ of the DNS SGS stress $\tau$ and the modelled SGS stress $\tau^{model}$ are used as a priori tests to evaluate the performance of the SGS models. $C(\tau)$ is defined as

where $\left\langle{~{}*~{}}\right\rangle$ denotes the averaging over the spatial volume.

Fig. 14 presents the RMSE and correlation coefficients of the GEP, SRNN, and GEPNN SGS models during training progress. For each method, we run $30$ training sessions, and the lines in Fig. 14(a) represents the best results for each method, while the shaded areas indicate variations among different training sessions. It is firstly noted that it is quite difficult for the SRNN2 with its simple structure to converge to a reasonable result, even using the same hyperparameters as for the SRNN1. This illustrates the limitation of the inflexible preset structure of SRNN. Furthermore, the RMSE and correlation coefficients obtained by the SRNN1 and GEPNN models are both very promising, better than the traditional GEP model, indicating that SRNN and GEPNN have better convergence performance given the same training time. Moreover, the GEPNN performs slightly better compared to the SRNN1, especially for the correlation coefficients shown in Fig. 14(b).

The final expressions of the GEP, SRNN1 and GEPNN models, which are selected as the best one from the $30$ training sessions for each method, are as follows:

It is noted that $g^{(1\dots 4)}$ contain no invariants $I_{1\dots 4}$, and all terms of $g^{(1\dots 4)}$ are constants. Table 4 summarises the constant coefficients $g^{(1\dots 4)}$ for the GEP, SRNN1 and GEPNN models obtained from $30$ training sessions. The optimization process of GEP includes a series of random processes, for example, the probability of mutation and crossover. Thus, the convergence results of the GEP method may differ for each training run within a limited time. As a comparison, the GEPNN shows good convergence within a limited training time because the neural network optimization step is very efficient with gradient information to obtain the correct optimization direction. Among the $30$ training sessions of SRNN1, however, the trained coefficients fucntions $g^{(1)}$ and $g^{(3)}$ are either $0$, or complex combinations of invariants $I_{1\dots 4}$ with coefficients of small magnitudes (not shown in Table 4). This suggests that the strategies we introduced to enhance the model sparsity in GEPNN do result in better convergence to simplified models compared to the original SRNN.

Furthermore, we apply the trained GEP and GEPNN models to LES of homogeneous isotropic turbulence at a grid resolution of $128^{3}$, and the a posteriori results from these ML models are then investigated. Moreover, the traditional SGS models, i. e. DSM and DMM, are also applied and shown for comparison. Fig. 15 shows the spectrum of the velocity field. The velocity spectrum of the DNS has a long inertial region. Owing to the filtering operation, the fDNS case deviates in the high-wavenumber region. For comparison, we conducted an LES without SGS models at a grid resolution of $128^{3}$, i.e. an under-resolved DNS. A significant deviation occurred near the cutoff wavenumber because the simulation failed to account for the dissipation of the unresolved scales. The spectra of the DSM and DMM models, however, are damped near the cutoff wavenumber but are energy-rich at lower wavenumbers. The GEP and GEPNN models predict velocity spectra that nearly overlap with the fDNS data, with the GEPNN model being more accurate near the cutoff wavenumber.

The SGS energy flux is also evaluated to investigate the energy transfer between the resolved and subgrid scales, which is defined as the product of the SGS stress tensor and the strain rate tensor

Fig. 16 shows the PDFs of the SGS energy flux $\Pi$ for the LES using different models, where the positive SGS flux refers to the energy transfer from the resolved field to the subgrid field, and the negative flux is the energy transfer in the opposite direction. Compared to the DSM and DMM models, the GEP model predicts a positive SGS flux similar to the fDNS data, but overpredicts the negative flux. However, the GEPNN model makes a precise prediction, not only for the forward energy transfer but also for the backscatter of kinetic energy from the subgrid field to the resolved scale.

The contours of one component of the SGS stress tensor $\tau_{ij}^{A}$ at an arbitrary slice of the fDNS and LES using different SGS models are shown in Fig. 17, where regions with high stress intensities represent localized subgrid stress structures. The structures generated by the DMM model are larger and more diffusive than the fDNS data and those generated by data-driven models, because the DMM introduces extra dissipation that dissipates small-scale structures. However, the small-scale structures can be effectively captured using the GEP and GEPNN models, similar to the fDNS data, illustrating that data-driven closures of SGS stress can model small-scale flow behaviour better than the DMM model. This observation is consistent with the previous findings and demonstrates the superior performance of the data-driven SGS models, as illustrated in Fig. 15 and Fig. 16.

The Reynolds number is a dimensionless variable defined as the ratio of inertial force to viscous force, and a high Reynolds number indicates that the influence of inertia is significant. As the flow properties vary across different Reynolds numbers, the generalization ability of the SGS model at different Reynolds numbers is critical. The spectrum of the velocity fields obtained from LES at different Reynolds numbers and grid resolutions are presented in Fig. 18. This shows that data-driven models with functional expressions have excellent generalization ability, allowing them to adapt to different flow conditions.

In the SGS stress modelling problem, the GEPNN showed better convergence performance than traditional GEP and SRNN methods. In particular, using a simpler structure, SRNN2 failed to model the SGS stress, and all the training sessions of SRNN2 diverged. This phenomenon reflects that the preset structure of SRNN has a strong influence on performance. Ultimately, a reliable structure needs more experience and analysis for the modelling problem, which becomes an obstacle for applying SRNN in unknown fields. The GEPNN model achieved the lowest RMSE and highest correlation coefficient in the a priori test, which shows that the predictions of the GEPNN model correlated more closely with the DNS data. In the a posteriori test, LES using the GEPNN model could predict turbulence statistics that align with the filtered DNS data. Compared to the traditional SGS models, the turbulence structures predicted using the GEPNN model are much finer. Moreover, the GEPNN model has a strong generalization ability, which means that the algebraic model trained from one flow condition can perform well for other Reynolds numbers and grid resolutions.