#  Turbulence closure modeling with data–driven techniques:
Investigation of generalizable deep neural networks

Salar Taghizadeh1<EMAIL_ADDRESS>Freddie D. Witherden2 Yassin A.
Hassan1,3 Sharath S. Girimaji1,2,4 1J. Mike Walker ’66 Department of
Mechanical Engineering, Texas A&M University, College Station, TX 77843
2Department of Ocean Engineering, Texas A&M University, College Station, TX
77843  3Department of Nuclear Engineering, Texas A&M University, College
Station, TX 77843  3Department of Aerospace Engineering, Texas A&M University,
College Station, TX 77843

###### Abstract

Generalizability of machine–learning (ML) based turbulence closures to
accurately predict unseen practical flows remains an important challenge. It
is well recognized that the neural network (NN) architecture and training
protocol profoundly influence the generalizability characteristics. At the
Reynolds–averaged Navier–Stokes (RANS) level, NN–based turbulence closure
modeling is rendered difficult due to two important reasons: inherent
complexity of the constitutive relation arising from flow–dependent
non–linearity and bifurcations; and, inordinate difficulty in obtaining
high–fidelity data covering the entire parameter space of interest. Thus, a
predictive turbulence model must be robust enough to perform reasonably
outside the domain of training. In this context, the objective of the work is
to investigate the approximation capabilities of standard moderate–sized
fully–connected NNs. We seek to systematically investigate the effects of:
_(i)_ intrinsic complexity of the solution manifold; _(ii)_ sampling procedure
(interpolation vs. extrapolation) and _(iii)_ optimization procedure. To
overcome the data acquisition challenges, three proxy–physics turbulence
surrogates of different degrees of complexity (yet significantly simpler than
turbulence physics) are employed to generate the parameter–to–solution maps.
Lacking a strong theoretical basis for finding the globally optimal NN
architecture and hyperparameters in the presence of non–linearity and
bifurcations, a ‘brute–force’ parameter–space sweep is performed to determine
a locally optimal solution. Even for this simple proxy–physics system, it is
demonstrated that feed–forward NNs require more degrees of freedom than the
original proxy–physics model to accurately approximate the true model even
when trained with data over the entire parameter space (interpolation).
Additionally, if deep fully–connected NNs are trained with data only from part
of the parameter space (extrapolation), their approximation capability reduces
considerably and it is not straightforward to find an optimal architecture.
Overall, the findings provide a realistic perspective on the utility of ML
turbulence closures for practical applications and identify areas for
improvement.

##  I Introduction

Turbulent flows exhibit vastly different characteristics in different
parameter regimes depending upon the mean strain–to–rotation rate ratio
blaisdell1996simulation ; girimaji2000pressure , mean-to-turbulence time scale
ratio speziale1991analytical ; mishra2010pressure ; mishra2017toward ,
underlying flow instabilities speziale1996consistency , large–scale
unsteadiness, presence or absence of system–rotation or streamline–curvature
speziale1991analytical , body–force effects and flow geometry. In addition,
due to inherent complexity of the turbulence phenomenon, closure models
developed in one parameter regime cannot be presumed to be reasonable or even
valid in other regimes. The strong dependence of flow statistics on the
various physical parameters is one of the enduring challenges in the field of
turbulence closure model development. The degree of difficulty of closure
modeling depends upon the level of closure. In the Reynolds–averaged
Navier–Stokes (RANS) method, the fundamental governing equations are averaged
over all scales of motion leading to significant reduction in computational
effort needed for performing a flow simulation. The reduction in computational
burden comes at the cost of increased complexity of closure modeling. Ad hoc
simplifications or assumptions are typically invoked to close various terms in
the RANS equations. While RANS models may perform adequately in flows in which
they are calibrated, they can be catastrophically wrong in other complex flows
due to lack of generalizability. Despite inherent limitations, RANS is widely
used in industrial applications involving complex flows due to ease of
computations. At the other extreme of the closure spectrum, in the large–eddy
simulation (LES) approach, dynamically important scales are resolved and only
the small–scale motions are modeled. The small–scales are significantly easier
to model as they embody most of the ‘universal’ aspects of turbulence and
therefore, are more easily amenable to generalizability than their RANS
counterpart. Thus, in LES, the relative simplicity and generalizability of
subgrid closure models comes at the expense of significantly increased
computational costs. The closure modeling challenges of scale resolving
simulations (SRS) are of intermediate degree of difficulty as they resolve
more scales than RANS but significantly lesser than LES girimaji2006partially
; girimaji2006partially2 .

There is heightened expectation in recent times that ML techniques can be used
to significantly improve RANS turbulence closure model performance in complex
industrial flows. The rationale is that the shortcomings of the physics–based
closures can be circumvented by appropriate data–based training of the models.
Toward this end, many authors have used different ML methods to model the
turbulence constitutive relation at all level of closure modeling, LES, SRS
and RANS. ML–enhanced LES closures have been proposed in different studies
sarghini2003neural ; gamahara2017searching ; maulik2017neural ; maulik2018data
; beck2019deep ; zhou2019subgrid ; xie2020modeling ;
schoepplein2018application ; reissmann2021application ; maulik2019subgrid ;
wang2018investigations ; maulik2019sub ; xie2019artificial ; xie2020artificial
; xie2019modeling ; park2021toward ; wang2021artificial ; subel2021data ;
pawar2020priori . ML–enhanced RANS closures have also been proposed in
numerous works using different ML algorithms including NNs ling2016reynolds ;
zhang2019application ; sotgiu2019towards ; geneva2019quantifying ;
fang2020neural ; jiang2020novel ; yin2020feature ; taghizadeh2020turbulence ;
jiang2021interpretable ; zhou2021learning ; yang2020improving ;
volpiani2021machine ; maulik2021turbulent ; zhu2019machine ; peters2021s ,
random forest breiman2001random ; wu2017priori ; kaandorp2020data ;
wu2018physics ; yang2020improving , gene expression programming (GEP)
weatheritt2017development ; weatheritt2017machine ; akolekar2018development ;
zhao2020rans ; lav2019framework ; waschkowski2021multi  and deterministic
symbolic regression models schmelzer2020discovery ; huijing2021data ;
beetham2020formulating ; zhang2021customized . A complete list of important
contributions in this area can be found in recent review papers kutz2017deep ;
duraisamy2019turbulence ; brunton2020machine ; beck2021perspective ;
duraisamy2021perspectives . ML methods for turbulence closure are not without
their own challenges and shortcomings. There is no clear guidance on the
optimal choice of architecture and hyperparameters or different elements of
the training procedure. For instance, in order to develop a new algebraic
Reynolds stress model using channel flow dataset, different NN architectures
have been employed in literature. Zhang et al. zhang2019application  trained a
NN with 4 hidden layers and 20 neurons per layer. Fang et al. fang2020neural
employed a NN with 5 hidden layers and 50 neurons per layer. Jiang et al.
jiang2020novel  used a DNN with 9 hidden layers and varying number of neurons
in each layer. For planner and periodic hill channel flow dataset, Sotgiu et
al. sotgiu2019towards  used a NN with 8 hidden layers and 8 neurons in each
layer. Lacking a formal procedure for NN selection and training protocol,
ML–assisted turbulence closures can be as ad hoc as the traditional models
and, more importantly, lack generalizability. In other areas, it has been
shown that ML models developed with subject matter expertise has a better
chance of succeeding anand2020black . Thus, it is of much interest to examine
generalizability in turbulence–like systems that are much simpler to compute.

The earliest fundamental analyses of the approximation capabilities of NNs
demonstrate that any continuous function on a bounded domain can be
approximated to arbitrary accuracy with at least one hidden layer with
sufficiently many neurons cybenko1989approximation ; hornik1989multilayer .
However, these results do not quantify the required sizes of NNs to achieve
these rates. Approximation rates of NNs with different activation functions
for various functions are presented in later studies barron1993universal ;
leshno1993multilayer ; maiorov1999lower ; mhaskar1993approximation ;
mhaskar1996neural ; yarotsky2017error ; petersen2018optimal . As the most
efficient statistical learning methods, deep neural networks (DNNs) have been
used to directly model the solution of PDEs raissi2019physics ; jin2021nsfnets
; berner2020analysis ; weinan2017deep ; han2018solving ; sirignano2018dgm  in
which it is often plausible to bound the size of the involved NNs in a way
that overcome the curse of dimensionality, i.e., the approximation rates do
not grow exponentially with increasing input dimension. DNNs have also been
used to approximate the parameter (solution) maps in parametric problems
dal2020data ; hesthaven2018non ; khoo2021solving ; lee2020model ;
kutyniok2021theoretical ; geist2021numerical . In parametric PDE problems,
parametric map connects the solution space of PDE to the parameters that
describe mathematical and physical constraints of PDE, for instance shape of
physical domain, boundary conditions. For linear parametric PDEs it has been
theoretically shown that feed–forward NNs of sufficient depth and size with
Rectified Linear Unit (ReLU) activation function are able to produce very
efficient approximations. However, the results are not optimal and do not
yield minimum depth and minimum number of neurons per layer, thus it is not
clear whether deep NNs are indeed necessary kutyniok2021theoretical .

The practical learning problems consists of several aspects: _(i)_ capacity of
NN architecture in describing the data; _(ii)_ data availability to describe
the true model and _(iii)_ optimization procedure in finding the best fit
kutyniok2021theoretical . In general, deep learning is a non–convex
optimization problem and while a good approximation of a given unknown
function by a NN may exist, it is unclear how that can be expeditiously
determined in a practical application. Additionally, it is certainly feasible
that the generated data might not contain sufficient information from the true
model to assure that the optimization process will converge to the
theoretically best approximation geist2021numerical . Therefore, it is
necessary to numerically examine the existence of generalizable,
reasonably–sized NNs from practical learning prospective. In this work we use
supervised learning method, to train standard fully–connected NNs to fit the
parameter–to–solution maps in the context of turbulence RANS closure modeling.

Turbulence has long been recognized as a complex flow phenomenon due to
inherent chaotic tendencies, multi–scale interactions arising out of the
non–linearity of the governing equations. As a consequence, a turbulence
constitutive relation which relates Reynolds stresses to the mean–flow field
statistics can exhibit multiple complex features:  _(i)_ strong history and
inhomogeneity effects;  _(ii)_ non–linearity of unknown degree and  _(iii)_
multiple bifurcations in complex flows involving compressibility, system
rotation, streamline curvature, additional body forces and complicated
geometric features. Further, the RANS closures must satisfy key conservation
laws, physical principles and mathematical constraints
taghizadeh2020turbulence . Thus, generating a reasonable parameter–to–solution
map for turbulence modeling at the RANS level is challenging. Even if the full
parameter space is identified, obtaining data from direct numerical
simulations (DNS) or experiments could be prohibitively expensive and indeed
infeasible.

At the current state of development, data–driven methods work best for
interpolation, while application to parameter space outside the training
domain is fraught with uncertainty. When a NN–based turbulence model is used
for computing flows of engineering interest, parameter regimes outside the
training domain will likely be encountered. Therefore, it will be useful to
assess the extrapolation capabilities of the NN–based models.

To circumvent data sampling limitations, we propose to employ simplified
proxy–physics turbulence surrogates to generate parameter–to–solution maps.
Using this data, we systematically examine three main challenges of practical
learning problems in the context of turbulence closure modeling:

  1. 1.

_Effect of the intrinsic complexity:_ The effect of the complexity of the
solution manifold on the performance of the NNs is investigated by employing
three proxy–physics surrogates of different degrees of non–linearity and
bifurcation characteristics. The turbulence surrogates mimic some of the key
features of turbulence and provide sufficiently many training/testing data at
low computational expense. It must be iterated that the actual turbulence
constitutive relationship can be considerably more complex due to history
(transient) effects, flow inhomogeneity and multiple bifurcations, which are
not considered in the proxy–physics models.

  2. 2.

_Effects of sampling procedure (interpolation vs. extrapolation):_ The
proxy–model approach allows for training data to be generated over all of the
parameter space or select regimes to investigate the effects of sampling.
Based on underlying physics of these proxy–physics models, training and
testing datasets are segmented into three different scenarios to assess
generalizability (interpolation vs. extrapolation) characteristics.

  3. 3.

_Effects of the optimization procedure:_ Architecture and the hyperparameters
of the DNNs can significantly affect their performance and generalizability
characteristics. In this study we examine the existence of generalizable NNs
by performing a systematic search in state–space of hyperparameters and
network–size for different data sampling scenarios. We seek to establish the
optimal choice of hyperparameters and DNN neurons of required to achieve a
predetermined level of accuracy.

To enable a reasonably rigorous analysis, we _(i)_ restrict our consideration
to statistically two–dimensional homogeneous turbulence, which represents the
most elementary non–trivial flows of interest; and _(ii)_ utilize data
generated from proxy–physics turbulence surrogates that incorporate some of
the key aspects of real homogeneous flows. We consider using standard
fully–connected NNs to reproduce the results of the simplified proxy–physics
models based around a cubic polynomial. Such a model represents a gross
oversimplification of the true physical dynamics of turbulence. Hence, if a
moderate–sized NN, given an arbitrarily large amount of training data is
unable to adequately reproduce this model, it would be unreasonable to expect
NN–based ML models to perform well in unseen turbulent flows. Further, the
lessons learned for this simpler system can yield valuable insight into the
closure modeling of the more complex turbulence phenomenon.

The various NN hyperparameters and training elements considered in this study
are: number of hidden layers or depth of network, number of neurons in each
layer or width of network, type of loss function, type of activation function
for neurons, type of training optimizer, learning rare, batch size and
regularization coefficient. The organization in the rest of the paper is as
follows. The selected proxy–physics turbulence surrogates are explained in
Sec. II. The NN hyperparameters and training elements are described in Sec.
III. Section IV details how data is generated with the proxy–physics
turbulence surrogates in different regimes of two–dimensional homogeneous
turbulence. The results from different architectures and training methods are
presented in Sec. V. We conclude with a summary of findings and
recommendations in Sec. VI.

##  II Proxy–Physics Methodology

In a statistically steady two–dimensional homogeneous incompressible turbulent
flow field without body forces, only two parameters govern the physical
behavior – normalized mean strain ($s_{ij}$) and rotation rates ($r_{ij}$),

$s_{ij}=\frac{k}{\epsilon}S_{ij},\quad r_{ij}=\frac{k}{\epsilon}R_{ij},$ (1)

where

$S_{ij}=\frac{1}{2}(\frac{\partial U_{i}}{\partial x_{j}}+\frac{\partial
U_{j}}{\partial x_{i}}),\quad R_{ij}=\frac{1}{2}(\frac{\partial
U_{i}}{\partial x_{j}}-\frac{\partial U_{j}}{\partial x_{i}}).$ (2)

Here $k$ and $\epsilon$ are respectively turbulent kinetic energy and
turbulent dissipation. This set of flows do not exhibit any large–scale
instabilities, coherent structures, complex geometrical features or
statistical unsteadiness. Despite the apparent simplicity, homogeneous
two–dimensional turbulent flows embody multiple complex phenomena that are
strongly dependent upon the parameter values. The flow goes from hyperbolic to
rectilinear to elliptic streamline geometry with increasing mean rotation
rate. Depending upon the mean flow to turbulence strain rate ratio, the flow
physics can range from the rapid distortion limit to decaying anisotropic
turbulence.

Using representation theory, the four–term expansion of the anisotropy tensor
in terms of the normalized strain rate ($s$) and normalized rotation rate
($r$) tensors for two–dimensional mean flows can be expressed as
girimaji1996fully ,

$b_{ij}=G_{1}(s_{ij})+G_{2}(s_{ik}r_{kj}-r_{ik}s_{kj})+G_{3}(s_{ik}s_{kj}-\frac{1}{3}\delta_{ij}s_{mn}s_{nm})+G_{4}(r_{ik}r_{kj}-\frac{1}{3}\delta_{ij}r_{mn}r_{nm}),$
(3)

where scalar coefficients $G_{1}-G_{4}$ are constitutive closure coefficients
(CCC) taghizadeh2020turbulence  that must be modeled. They are functions of
scalar invariant of the strain and rotation–rate tensors
($\eta_{1}=s_{ij}s_{ij},\eta_{2}=r_{ij}r_{ij}$) and other flow quantities such
as - turbulent kinetic energy ($k$), turbulence frequency ($\omega$) and the
coefficients of the pressure–strain correlation model.

As mentioned in the Introduction, generalizability of the NN models in this
simple system is a prerequisite to generalizability in actual turbulent flows.
Further, due to the fact that proxy–physics turbulence surrogates capture many
key statistical features of turbulence, this study will yield much valuable
insight into RANS ML modeling. Complete knowledge of the proxy–physics
solution enables precise error assessment incurred during generalization to
test flows. In this study, the proxy–physics data in different parts of the
domain are used to train NNs of different architectures and hyperparameters.
The networks are then tested in other parameter regimes to examine
generalizability.

###  II.1  Proxy–Physics Surrogates

As stated before the main objective of this work is to look at the challenges
posed by non–linearity and bifurcation effects of turbulence, therefore three
proxy–physics surrogates of different non–linearity and bifurcation
characteristics are employed to generate the parameter–to–solution maps: two
algebraic Reynolds stress models (ARSM) and one non–linear constitutive
relationship.

_Algebraic Reynolds Stress Model._ ARSM does reasonably well in capturing key
flow physics in different parts of the flow regime gatski1993explicit ;
girimaji1996fully . The model requires the solution of a cubic equation and
appropriate root must be chosen in different flow regimes to yield the correct
behavior girimaji1996fully . In this study, we use three–term self–consistent,
nonsingular and fully explicit algebraic Reynolds stress model (EARSM)
proposed by Girimaji girimaji1996fully  with two pressure–strain correlation
models, Launder, Reece and Rodi (LRR) launder1975progress  and and Speziale,
Sarkar and Gatski (SSG) speziale1991modelling  to generate stress–strain
(constitutive) relationship datasets for different normalized strain and
rotation rates.

An implicit algebraic equation for the anisotropy tensor can be obtained in
the weak equilibrium limit of turbulence using the following simplification
girimaji1996fully ,

$\frac{Db_{ij}}{Dt}=\frac{\partial b_{ij}}{\partial t}+U_{k}\frac{\partial
b_{ij}}{\partial x_{k}}\approx 0,$ (4)

where $D/Dt$ is the substantial derivative following the mean flow. The
weak–equilibrium assumption is valid for many flows wherein the timescale of
anisotropy evolution is rapid compared to the timescales of mean flow,
turbulent kinetic energy, and dissipation rate girimaji2001lower ;
gomez2014explicit . Using Eqs. (3) and (4), the non–linear algebraic Reynolds
stress equation with the three–term model ($G_{1}-G_{3}$) can be written as
the following cubic fixed–point equation for $G_{1}$ girimaji1996fully ,

$(\eta_{1}L_{1}^{1})^{2}G_{1}^{3}-(2\eta_{1}L_{1}^{0}L_{1}^{1})G_{1}^{2}+\Big{[}(L_{1}^{0})^{2}+\eta_{1}L_{1}^{1}L_{2}-\frac{2}{3}\eta_{1}(L_{3})^{2}+2\eta_{2}(L_{4})^{2}\Big{]}G_{1}-L_{1}^{0}L_{2}=0.$
(5)

Then $G_{2}$ and $G_{3}$ can be expressed as girimaji1996fully ,

$G_{2}=\frac{-L_{4}G_{1}}{L_{1}^{0}-\eta_{1}L_{1}^{1}G_{1}},\quad
G_{3}=\frac{2L_{3}G_{1}}{L_{1}^{0}-\eta_{1}L_{1}^{1}G_{1}}.$ (6)

The closure coefficients in Eqs. (5) and (6) are as follows:

$L_{1}^{0}=\frac{C_{1}^{0}}{2}-1,\quad L_{1}^{1}=C_{1}^{1}+2,\quad
L_{2}=\frac{C_{2}}{2}-\frac{2}{3},\quad L_{3}=\frac{C_{3}}{2}-1,\quad
L_{4}=\frac{C_{4}}{2}-1,$ (7)

where the $C$’s are numerical constants of the pressure–strain correlation
models. The numerical constants for LRR pressure–strain correlation model,
linear in the anisotropy tensor ($C_{1}^{1}=0$), and SSG, quasilinear in the
anisotropy tensor ($C_{1}^{1}\neq 0$), are given in Table 1.

Table 1: Coefficients in the LRR and SSG models Model |  $C_{1}^{0}$ |  $C_{1}^{1}$ |  $C_{2}$ |  $C_{3}$ |  $C_{4}$
---|---|---|---|---|---
LRR |  3.0 |  0 |  0.8 |  1.75 |  1.31
SSG |  3.4 |  1.8 |  0.36 |  1.25 |  0.4

The cubic relation in Eq. (5) has multiple real and complex roots and the
selection of the appropriate solution for $G_{1}$ is not straightforward. By
considering two physical selection criteria; _(i)_ continuity of $G_{1}$; and
_(ii)_ $G_{1}\leq 0$, Girimaji girimaji1996fully  derived a fully explicit
solution of the cubic Eq. (5) for scalar coefficient $G_{1}$ as follows:

$G_{1}=\left\\{\begin{array}[]{ll}\frac{L_{1}^{0}L_{2}}{(L_{1}^{0})^{2}+2\eta_{2}(L_{4})^{2}}&\mbox{for
}\eta_{1}=0,\\\
\frac{L_{1}^{0}L_{2}}{(L_{1}^{0})^{2}-\frac{2}{3}\eta_{1}(L_{3})^{2}+2\eta_{2}(L_{4})^{2}}&\mbox{for
}L_{1}^{1}=0,\\\
-\frac{p}{3}+(-\frac{b}{2}+\sqrt{D})^{\frac{1}{3}}+(-\frac{b}{2}-\sqrt{D})^{\frac{1}{3}}&\mbox{for
}D>0,\\\ -\frac{p}{3}+2\sqrt{\frac{-a}{3}}cos(\frac{\theta}{3})&\mbox{for
}D<0,b<0,\\\
-\frac{p}{3}+2\sqrt{\frac{-a}{3}}cos(\frac{\theta}{3}+\frac{2\pi}{3})&\mbox{for
}D<0,b>0.\end{array}\right.$ (8)

Here the discriminant $D$ of the cubic Eq. (5) is calculated as:

$D=\frac{b^{2}}{4}+\frac{a^{3}}{27},$ (9)

other parameters of Eqs. (8) and (9) are defined as below:

$\begin{split}a=(q-\frac{p^{2}}{3}),&\quad b=\frac{1}{27}(2p^{3}-9pq+27r),\\\
p=-\frac{2L_{1}^{0}}{\eta_{1}L_{1}^{1}},\quad
q=\frac{1}{(\eta_{1}L_{1}^{1})^{2}}&\Big{[}(L_{1}^{0})^{2}+\eta_{1}L_{1}^{1}L_{2}-\frac{2}{3}\eta_{1}(L_{3})^{2}+2\eta_{2}(L_{4})^{2}\Big{]},\\\
r=-\frac{L_{1}^{0}L_{2}}{(\eta_{1}L_{1}^{1})^{2}},&\quad
cos(\theta)=\frac{-b/2}{\sqrt{-a^{3}/27}}.\end{split}$ (10)

The first two cases in Eq. (8) are special limiting cases of the last three
and it can be shown that the limiting behavior can be easily calculated from
the general expressions girimaji1996fully . After the coefficient $G_{1}$ is
calculated, other CCC can also be obtained using Eq. (6) in the entire
parameter space.

The goal of this study is to use proxy–physics turbulence surrogates to
provide some turbulence subject matter expertise. For this reason, it is
necessary to ensure that the selected proxy model adequately incorporates some
of the known features of turbulence. Different turbulence states covered by
the ARSM turbulence surrogate are depicted for the considered parameter space
in Fig. 1 ($S\equiv\sqrt{s_{ij}s_{ij}}\equiv\sqrt{\eta_{1}}$ and
$R\equiv\sqrt{r_{ij}r_{ij}}\equiv\sqrt{\eta_{2}}$). At rapidly strained
turbulence state, strain rate dominates over rotation rate and the
discriminant $D$ (Eq. (9)) of the cubic fixed–point equation for $G_{1}$ is
negative. The realizability violations occur at this rapidly strained region
due to the governing elastic constitutive relationship. Rotation rate
dominates over strain rate at high values of mean vorticity state. Other
important turbulence states are also represented in this figure. Hyperbolic
streamline flows occur when $S>>R$; and the streamlines are elliptic when
$S<<R$. When $S\approx R$ rectilinear shear flows are seen
mishra2013intercomponent .

Figure 1:  Different states of turbulence of ARSM

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(g)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(h)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(i)

Figure 2:  Parameter–to–solution maps, (a), (d) and (g) ARSM with SSG, (b),
(e) and (h) ARSM with LRR, (c), (f) and (i) SZL

The parameter–to–solution maps obtained from ARSM with two pressure–strain
correlation models are shown in Fig. 2. It is shown that the coefficient
$G_{1}$ is a continuous function across bifurcation line $D=0$ in ARSM model.
Therefore, $G_{1}$, i.e., the effective turbulent viscosity is well defined in
the entire parameter space wherein different important turbulence states are
covered. Fig. 2 shows well defined values for $G_{2}$ and $G_{3}$ coefficients
in the entire parameter space. The corresponding CCC contours for ARSM with
both the pressure–strain correlation models have similar shapes. However, CCC
values are in wider ranges with SSG model. In particular, the magnitude of the
$G_{1}$ is approximately zero, $10^{-5}$ for larger values of parameters,
while it is in the order of $10^{-1}$ in decaying turbulence state. As
mentioned earlier, small values of $G_{1}$ for large $S$ is a consequence of
elastic constitutive behaviour and is needed for preserving realizability.

_Non–linear Constitutive Relationship._ In addition to the three–term ARSM
surrogate proposed by Girimaji girimaji1996fully  with different
pressure–strain correlation models, the four–term non–linear constitutive
relationship proposed by Shih et al. shih1993realizable (Shih, Zhu and Lumley
– SZL model) is also considered as a proxy–physics turbulence surrogate to
generate the parameter–to–solution map in this study. In SZL model, the CCC
are expressed as shih1993realizable ; craft1996development ,

$\begin{split}G_{1}&=-C_{\mu}=\frac{-\frac{2}{3}}{1.25+\sqrt{2\eta_{1}}+0.9\sqrt{2\eta_{2}}},\quad
G_{2}=\frac{-\frac{15}{2}}{1000+(\sqrt{2\eta_{1}})^{3}},\\\
G_{3}&=\frac{\frac{3}{2}}{1000+(\sqrt{2\eta_{1}})^{3}},\quad\hskip
78.24507ptG_{4}=\frac{\frac{-19}{2}}{1000+(\sqrt{2\eta_{1}})^{3}}.\end{split}$
(11)

It should be noted that the CCC expressions proposed in SZL model are simpler
compared to ARSM model counterparts. The non–linear SZL model does not have
bifurcation and the continuous function of $G_{1}$ is obtained with one single
relation as given by Eq. (11) in the entire parameter space. However, in ARSM
two different relationships are employed. For the considered parameter space,
the third relation in Eq. (8) is used in the rotation–dominated region $D>0$
and the fifth relation is used in the strain–dominated region $D<0$.
Therefore, the parameter–to–solution map generated by SZL model has less
complexity level compared to ARSM, Fig. 2. In SZL model, $G_{2}-G_{4}$
coefficients are only functions of the invariant of strain rate ($\eta_{1}$),
therefore contour plots show vertical lines in the parameter space. Overall,
the physics underlying CCC contours is a reasonable facsimile of real
turbulent flows at different complexity levels of proxy–physics turbulence
surrogates.

Therefore, three proxy–physics turbulence surrogates of different degrees of
complexity can be used to examine the challenges posed by non–linearity and
bifurcation effects of the solution manifold: _(i)_ SZL model, non–linear
constitutive relation with no bifurcation in the parameter space; _(ii)_ ARSM
with LRR model, mildly non–linear constitutive relation with bifurcation in
the regime of interest; and _(iii)_ ARSM with SSG model, moderately non–linear
constitutive model with bifurcation in the parameter space. The objective of
this study is to investigate if a reasonably–sized, fully–connected NN, given
an arbitrarily large amount of training data can simulate the simple
surrogates of turbulence. If the NN–based ML models can not perform well with
the simplified descriptions, it is unlikely to perform successfully in real
unseen turbulent flows.

##  III  Machine Learning

_Previous investigations._ ML is a general term to describe a class of
algorithms which uses data to generate models. The selection of modeling
strategy depends on the type of problem. Recently, DNNs have been widely used
for modeling turbulent flows. Ling et al. ling2016reynolds  trained a DNN with
8 hidden layers, and 30 neurons per hidden layer with available high–fidelity
data: duct and channel flow, perpendicular and inclined jet in cross–flow,
flow around a square cylinder and flow through a converging–diverging flow.
Their results showed that incorporating the Galilean invariance property in
the DNN architecture can improve the predictive capability of ML turbulence
models. Zhang et al. zhang2019application  trained a NN with 4 hidden layers
and 20 neurons per layer in order to develop a model to predict the Reynolds
stress of a channel flow at different Reynolds numbers. They obtained well
behaved models by introducing regularization in their training algorithm.
Using same channel flow datasets, Fang et al. fang2020neural  trained a NN
with 5 hidden layers and 50 neurons per layer in order to develop an improved
Reynolds stress tensor model. Jiang et al. jiang2020novel  used DNS of channel
flow at different Reynolds number in order to developed a new algebraic
Reynolds stress model by training a DNN with 9 hidden layers and varying
number of neurons in each layer as, 12, 18, 21, 27, 32, 35, 30, 28, and 27,
respectively. Sotgiu et al. sotgiu2019towards  developed a Reynolds stress
constitutive model by training a NN with 8 hidden layers and 8 neurons. They
used planner channel and periodic hill channel flow datasets for training the
ML algorithm. Geneva and Zabaras geneva2019quantifying  trained a NN with 5
hidden layers and tapered the number of neurons in the last two hidden layers
to prevent weights from being too small and improve training performance. The
required training data was generated by performing LES simulations of
different flows: converging–diverging channel, periodic hills, square duct,
square cylinder and tandem cylinders. The importance of judicious choice of
network architecture has been recognized in the field of turbulence modeling.
However, there have been no studies in literature to examine if the inherent
complexity of turbulence, sampling and training procedure pose any further
challenges.

_Objective._ Determining a suitable network architecture and training
hyperparameters of a DNN is an empirical task and it has been shown that there
is a strong positive correlation between the final performance of the trained
NNs and experience of the user in optimizing the hyperparameters
anand2020black . There is no rigorous guidance on the right choice of
architecture and different elements of the training procedure in developing
ML–assisted turbulence models. The objective of this study is to examine the
existence of generalizable, reasonably–sized, fully–connected NNs under full
and partial availability of the training datasets. Adopted proxy–physics
turbulence surrogates are used to easily generate the training data for all
flows in the parameter space. Based on underlying physics of the simplified
proxy–physics models, training and testing datasets are segmented at three
different scenarios to assess the effect of sampling procedure on optimization
capabilities and generalizability characteristics (interpolation vs.
extrapolation) of the DNNs of different architectures and hyperparameters.

Figure 3:  Profiles of different activation functions

_Deep neural networks._ A DNN is a class of model which transforms the input
parameters (features) through several layers of units (neurons). Each unit
(neuron) is connected by affine linear maps between units in successive layers
and then with non–linear (scalar) activation functions within units. Different
activation functions such as, rectified linear unit (ReLU), leaky ReLU,
exponential linear unit (ELU), Sigmoid and hyperbolic tangent
goodfellow2016deep  are shown in Fig. 3. The ReLU function which is the
popular choice in the ML literature is defined as:

$\sigma(z)=max(z,0),$ (12)

where, $z$ is the input to a neuron. In this study, we examine the effects of
employing different activation functions on training and performance of the
NNs.

It should be noted that DNNs consist of simple functions and their efficiency
emanates from the interactions between large number of hidden layers
goodfellow2016deep . Due to their flexible architecture and superior
performance in modeling non–linear and complicated relationships with
high–dimensional data, they have become a popular subset of ML approaches.
Although a variety of network structures, such as convolutional neural
networks (CNN), recurrent neural networks (RNN), or long short–term memory
(LSTM) networks have been proposed in the ML literature goodfellow2016deep ,
we will restrict ourselves to fully–connected architectures. Fig. 4 shows the
schematic of a deep feed–forward NN (also termed as a multi–layer perceptron
(MLP)) which is trained using back–propagation with gradient descent method.
The input layer, the hidden layers and the output layer are also shown in this
figure.

Figure 4:  Illustration of a deep feed–forward fully–connected NN

_Loss function._ The DNNs are trained by minimizing a loss (objective)
function, which measures the difference between the predicted output of the
model and labels (ground truth data). It has been shown that the success of ML
models depends on the formulation of the loss function used for optimization
of the model coefficients (weights and biases of neurons) during the ML
training process taghizadeh2020turbulence . The type of loss function is also
problem specific and need to be selected properly. The root mean squared error
(RMSE) is the commonly used loss function:

$RMSE=\sqrt{\frac{\sum_{i=1}^{N}(y_{i}-\hat{y}_{i})^{2}}{N}},$ (13)

where $\hat{y}_{i}$ denotes the ML prediction and $y_{i}$ is the true labeled
data; $N$ is the number of training data. Alternatively, the mean absolute
percentage error (MAPE) is also considered,

$MAPE=\frac{1}{N}\sum_{i=1}^{N}{\frac{\lvert y_{i}-\hat{y}_{i}\rvert}{\lvert
y_{i}\rvert}},$ (14)

The RMSE loss works well in most cases, while the MAPE loss is better for the
case where the output has a large range of function values cai2021deepm . In
order to find the appropriate type of loss function, the models trained with
MAPE and RMSE are compared in this study.

Adding a regularization term in the loss function formulation is one of the
common ways of avoiding overfitting. Overfitting occurs when the expressivity
of the ML model is too large for the complexity of the function it
approximates. Introducing a regularization term in the loss function
formulation shrinks the model coefficients towards zero, decrease the
complexity of the model and hence significantly reduces the variance
james2013introduction . The $L_{1}$ norm (Lasso Regression) and $L_{2}$ norm
(Ridge Regression) are the two common regularization methods
james2013introduction ,

$\begin{split}MAPE(y,y_{i})+\lambda\sum_{i=1}^{N}{\lvert\omega_{i}\rvert},\quad
L_{1}-norm\\\ MAPE(y,y_{i})+\lambda\sum_{i=1}^{N}{\omega_{i}^{2}},\quad
L_{2}-norm\end{split}$ (15)

where $\omega$ denotes the weight of each neuron and $\lambda$ is a positive
hyperparameter to determine the strength of regularization. For this work, to
control the overfitting in some experiments, first a comparison between two
regularization methods, $L_{1}$ and $L_{2}$ norms are made and then best
performing method is applied in all the hidden layers during the training of
the NNs.

_Network hyperparameters._ Architecture and all training hyperparameters of
the NNs need to be suitably specified in order to build a robust and effective
model that can generalize to unseen datasets. For instance, hyperparameters
include size of the networks (number of layers or depth, width or breadth of
each layer), formulation and type of the loss function, type of regularization
and optimization algorithm, type of activation function (act), batch size
(bs), learning rate (lr), type of initialization and etc. Hyperparameter
optimization is an empirical task and grid search is usually adopted, i.e.,
many networks with several different combinations of interval values of each
hyperparameter are trained and compared based on their accuracy and
generalization ability.

In principle, all of the hyperparameters in the ML algorithm can be varied and
they might have significant impact on the model performance. However, the
predictive capability and generalization of a NN is mostly controlled by its
architecture; the depth and breadth of network. In this study we perform a
systematic search in state–space of hyperparameters and network–size to train
NNs that are efficient and ensure a low generalization error in interpolation
and extrapolation cases. We consider a matrix of network sizes by varying the
depth for four values of 1, 3, 5, 7 and the width for four values of 3, 5, 7,
15. Other hyperparameters examined in this study are shown in Table 2.

Table 2: Different values of the considered hyperparameters Hyperparameter |  Value
---|---
Activation function (act) |  ReLU, ELU, leaky ReLU, Tanh, Sigmoid
Learning rate (lr) |  $1\times 10^{-6}$, $1\times 10^{-5}$, $1\times 10^{-4}$, $1\times 10^{-3}$, $1\times 10^{-2}$
Batch size (bs) |  25, 50, 100, 1000
Optimization algorithm (opt) |  Adam kingma2014adam , RMSProp
Regularization coefficient $(\lambda)$ |  0.01, 0.1, 0.2
Initialization function |  Xavier normal glorot2010understanding

##  IV  Data Generation with proxy–physics turbulence surrogates

The data needed for training the ML algorithm is generated using the
proxy–physics surrogates as discussed in Sec. II. The number of data points
that are non–uniformly extracted for each input parameter ($\eta_{1}$,
$\eta_{2}$) in the range [$10^{-1},10^{4}$] is 150. Therefore, the total
number of data points is 22,500 in the entire parameter space. In this work,
different investigations are conducted in order to address the challenges in
finding optimum NN architecture and hyperparameters when training is performed
as follows:

  1. 1.

Fully available data in the entire parameter space from ARSM and SZL models

  2. 2.

Partially available data only in one part of parameter space from ARSM model

  3. 3.

Partially available data only in a narrow band of parameter space from ARSM
model

These three investigations lead to important inferences about existence of the
moderate–sized NN approximation solutions for data manifolds with different
complexity levels and their generalizability (interpolation vs. extrapolation)
characteristics. We generate the training datasets in three different
scenarios as follows.

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

Figure 5:  Training (green circles) and testing (red circles) datasets for (a)
Case–1, (b) Case–2, (c) Case–3

###  IV.1  Training data over the entire parameter space

This case represents the ideal scenario in which the training data fully
describe the true model within the parameter space. Thus there is no physics
that is unseen by the ML model which is trained in supervised manner.
Sufficient data is gathered over strain, shear and rotational flow regimes. It
should be noted that providing the real turbulence data in the entire
parameter space, requires expensive DNS over a wide range of flows. In this
case, generalizability is expected to be trivially straightforward. The
optimum choice of hyperparameters and necessary number of neurons of DNNs
needed to have a sufficiently accurate approximation are investigated. As one
of the objectives of this study is to examine the effects of underlying
physics of the solution manifold (non–linearity and bifurcation effects) on
training and performance of the ML models, three parameter–to–solution maps
generated by proxy–physics turbulence surrogates discussed in Sec. II are
considered for this ideal case. For all the experiments in this case, the
generated data points in the parameters space are randomly split into 75% for
training (Training Data) and 25% for validation and testing the model (Testing
Data). The testing data is used for final evaluation of the models. Fig. 5(a)
represents randomly distributed data in the entire parameter space for this
case.

###  IV.2  Training data only in the strain–dominated region (D ¡ 0)

In this case, the training data is partially available and restricted to a
subset of the parameter space. Therefore, part of the physics in the data
manifold is seen by the NNs. Data points in the strain–dominated region $D<0$
(data generated with the fifth relation in Eq. (8)) are used for training and
data points in the rotation–dominated region $D>0$ (data generated with third
relation in Eq. (8)) are used for testing. Fig. 5(b) illustrates the segmented
training and testing datasets for this case. This represents an important
generalizability challenge (extrapolation) as training and testing are
performed in parameter regimes of distinctly different turbulence physics.
Existence of generalizable DNNs for this case is investigated for
parameter–to–solution maps generated by ARSM with SSG and LRR pressure–strain
correlation models.

###  IV.3  Limited training data in the shear–dominated region

In many instances training data is partially available only in a very narrow
region of the parameter space. In this regard, we examine the existence of a
generalizble NN trained with the limited dataset that covers shear–driven
turbulence physics. In this case, the parameter–to–solution map generated by
ARSM with SSG proxy–physics model is considered. As it can be seen from ARSM
equation, Eq. (8), shear flow represents the bifurcation region between strain
and rotation flows. An arbitrary zone near the bifurcation line $D=0$ is
defined as, $-1000<(\eta_{1}-0.7\eta_{2})<1000$. The data points inside this
region are used for training and data points outside this region are used for
testing. The segmented regions are represented in Fig. 5(c).

Table 3: Selected hyperparameters Hyperparameter |  Value
---|---
Learning rate (lr) |  $1\times 10^{-3}\sim 1\times 10^{-6}$
Batch size (bs) |  50
Optimization algorithm (opt) |  Adam
Initialization function |  Xavier normal glorot2010understanding

##  V  Results

###  V.1  RMSE loss function vs. MAPE loss function

It has been shown that the appropriate formulation of loss function and
optimal choice of the flow statistics contributing to the loss function impact
the success of the ML trained turbulence models taghizadeh2020turbulence . As
mentioned in Sec. III we examine the performance of the ML models trained with
different loss functions in this study. For this analysis the
parameter–to–solution map generated by ARSM with SSG proxy–physics model is
considered. First, the randomly generated dataset in the entire parameter
space is segmented as 75% for training and 25% for testing of the models.
Then, we train two DNNs with RMSE and MAPE loss functions without any
regularization. The selected fixed network architecture for both of the cases
has 7 hidden layers with 7 computation neurons in each layer. In this case,
the type of activation function is ReLU and all other hyperparameters of the
models are as shown in Table 3. As shown in the table, a variable learning
rate is employed for training the models.

The performance of the models trained with different loss functions are
reported with both MAPE and RMSE error metrics on training and testing
datasets in each column of Table 4. It can be seen that the network with MAPE
loss function has smaller training and testing errors compared to the network
with RMSE loss function. As mentioned in Sec. II, the CCC have a large range
of values over the parameter space. When the RMSE is used as the loss
function, the training and back–propagation process are mostly dominated by
the CCC with larger magnitudes, as the small value coefficients contribute
very little to the NN optimization process. Local error contours of the models
trained with different loss functions are shown in Figs. 6 and 7. It can be
seen that the ML algorithm trained with MAPE has smaller local errors in
different turbulence physics regions in the entire parameter space.

Comparing different error metrics, MAPE vs. RMSE for final performance of the
models in Table 4 and absolute error vs. MAPE for local error contours in
Figs. 6 and 7, clearly shows that MAPE metric has easily interpretable
presentation of final model performance in this study. Therefore, MAPE is
selected as the appropriate type of loss function and evaluation metric for
the remainder of the analysis.

Table 4: Performance of models trained with different loss functions | Errors with various evaluation metrics
---|---
ML model |  RMSE–training |  RMSE–testing |  MAPE–training |  MAPE–testing
RMSE–loss func. |  $9.28\times 10^{-3}$ |  $9.34\times 10^{-3}$ |  4.26 |  4.24
MAPE–loss func. |  $4.03\times 10^{-3}$ |  $3.90\times 10^{-3}$ |  0.043 |  0.044

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

Figure 6:  Absolute error contours of CCC for ML models trained with different
loss functions, (a)-(c) RMSE, (d)-(f) MAPE

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

Figure 7:  MAPE contours of CCC for ML models trained with different loss
functions, (a)-(c) RMSE, (d)-(f) MAPE

###  V.2  Case–1: Training data over the entire parameter space

In this experiment, three parameter–to–solution maps generated by SZL, ARSM
with SSG and LRR proxy–physics turbulence surrogates are considered. The
datasets are randomly distributed for training the NNs, i.e., 75% for training
and 25% for testing the models. Using the generated datasets, networks with
different architectures are trained. The type of activation function for all
computation neurons is ReLU and all other hyperparameters are fixed as shown
in Table 3 for this case. Fig. 8 demonstrates the training and testing MAPE
errors of all the 16 network architectures trained with different data
manifolds. In this figure number of hidden layers and number of neurons in
each layer are shown with horizontal and vertical axes, respectively. Although
for this case, $L_{2}$ regularization is not used during ML training, all the
16 networks have similar performance on both the training and testing datasets
and models test quite well without over/under fitting. It can be seen that
shallow NNs (networks with one layer) and DNNs with small width (networks with
three neurons in each layer) have the worst training and testing errors for
all the data manifolds in the interpolation case. Additionally, as the number
of the layers and neurons in each layer increase, the training and testing
errors decrease and for this case NN with largest degrees of freedom (7L–15N)
has the lowest training and testing errors. This is a surprising observation
that the NN-based models require very large networks (more degrees of freedom
than true proxy–physics model) to reduce errors to a reasonable level even for
this interpolation case. Therefore, DNNs are not efficient approximations for
the non–linear solution manifolds created by simplified proxy–physics
turbulence surrogates.

The approximation capabilities of the large NNs with the data manifolds of
different complexity levels are comparable in this case. However, for the
networks with reasonable degrees of freedom (total number of neurons less than
25), the errors are smallest in data manifold with no bifurcation (SZL model).
These networks have bigger errors in data manifold with moderate non–linearity
(ARSM with SSG model) compared to data manifold with mild non–linearity (ARSM
with LRR model). Local MAPE contours in the entire parameter space for the
networks with 7 hidden layers and different neurons trained with data manifold
generated by ARSM and SSG model are shown in Fig. 9. It is evident that by
increasing the number of neurons in each layer, error decreases in the entire
parameter space and for the network with large architecture (7L–15N) the
maximum error occurs mainly near the bifurcation line $D=0$ where the physics
of turbulence undergoes rapid change. It should be noted that ReLU activation
function is considered for all the NNs so far.

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

Figure 8:  Training and testing MAPE for Case–1, (a) and (d) SSG, (b) and (e)
LRR, (c) and (f) SZL

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(g)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(h)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(i)

Figure 9:  MAPE contours for Case–1, (a)-(c) 7L–5N, (d)-(f) 7L–7N, (g)-(i)
7L–15N

For Case–1, the effects of using ReLU and Sigmoid activation functions on
approximation capability of the NNs with different architectures are further
investigated. In this experiment all the hyperparameters are kept fixed as
shown in Table 3 and the data manifold generated by ARSM with SSG model is
considered. The training and testing MAPE errors of all the 16 network
architectures for both activation functions are illustrated in Fig. 10. It is
evident that the performance of the NNs with the considered activation
functions are significantly different. For both activation functions, networks
with small width have large training and testing errors. But small width
networks with Sigmoid activation function have smaller training and testing
errors. Unlike the networks with ReLU activation function, the training and
testing errors does not decrease as the number of the layers increases in
networks with Sigmoid activation function and these networks outperform in
shallow with large width architectures. It should be noted that Sigmoid
activation function involves expensive operations (exponentials, etc.)
compared to ReLU which is simply thresholded at zero (Eq. (12)). As
illustrated in Fig. 3, Sigmoid activation function takes an input ($z$) and
outputs a value ($\sigma$) in the range between 0 and 1. It is known that when
a neuron with Sigmoid activation function saturates at either tail of 0 and 1,
it ‘kills’ the gradient and the information is not transferred through the
neuron goodfellow2016deep . Hence the saturation of the neurons with Sigmoid
activation function could negatively affect the learning of the large
networks. On the other hand, a large gradients for ReLU neurons in shallow and
small width networks could unfavorably update the weights during the
back–propagation process. Therefore, the ReLU neurons can irreversibly ‘die’
during training since they can get knocked off the data manifold
goodfellow2016deep .

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

Figure 10:  Training and testing MAPE for Case–1, (a) and (c) ReLU, (b) and
(d) Sigmoid

###  V.3  Case–2: Training data only in the strain–dominated region (D ¡ 0)

In this experiment the existence of generalizable NNs is examined when
training data is partially available only in one side of parameter space
(extrapolation). The data points in the strain–dominated region $D<0$ are used
for training and the data points in the rotation–dominated region $D>0$ are
used for testing of the models. For this analysis, ReLU activation function
and $L_{2}$ norm regularization with $\lambda=0.1$ are used. All other
hyperparameters are fixed as shown in Table 3. A justification for the choice
of hyperparameters for this case is given in appendix, Sec. APPENDIX. Fig. 11
shows the training and testing errors of the 16 network architectures trained
with partially available data generated by ARSM with SSG and LRR models.
Unlike the interpolation case (Case–1), the ML models trained with limited
(biased) data show overfitting for large networks. It is seen from Fig. 11
that for both the data manifolds, the NNs with one and three hidden layers
have the worst performance in training dataset. As the number of hidden layers
increases to five, the capability of the NNs in capturing the non–linear
relationship between the input parameters and CCC labels increases. Although,
the training error reduces with increasing number of hidden layers from three
to five, the testing error increases. By further increasing the number of
hidden layers and neurons the performance of the model in training and testing
datasets oscillates. Comparing the performance of the networks trained with
different data manifolds in this case, shows that similar to the interpolation
case, NNs trained with moderately non–linear data manifold (ARSM with SSG
model) have bigger errors than NNs trained with mildly non–linear data
manifold (ARSM with LRR model). However, unlike the interpolation case,
finding a generalizable NN with a reasonable size is not straightforward when
training data partially describe the true proxy–physics surrogate in the
parameter space (only the strain–dominated region in this case). Local MAPE
contours of DNNs with 7 hidden layers and different neurons trained with
moderately non–linear data manifold are compared in Fig. 12. It is seen that
DNNs with large number of neurons has relatively smaller errors in both
training and testing regions of parameter space for all the CCC.

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

Figure 11:  Training and testing MAPE for Case–2, (a) and (c) SSG, (b) and (d)
LRR

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(g)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(h)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(i)

Figure 12:  MAPE contours for Case–2, (a)-(c) 7L–5N, (d)-(f) 7L–7N, (g)-(i)
7L–15N

###  V.4  Case–3: Limited training data in the shear–dominated region

In this scenario the training dataset covers a limited but important region of
the parameter space near the bifurcation line $D=0$. The data points in rest
of the parameters space are used for testing the ML models. For this case, the
data manifold is generated by ARSM with SSG model. The ReLU activation
function is employed for all computation neurons and all other hyperparameters
are fixed as shown in Table 3. The performance of NNs with 16 architectures in
the training and testing datasets are compared in Fig. 13. Similar to Case–2,
a $L_{2}$ norm regularization with $\lambda=0.1$ is used during training to
reduce the overfitting of the ML models in this case. Although, the training
error reduces by increasing number of hidden layers from three to seven in
networks with large width, the testing error oscillates. Similar to Case–2,
selecting a reasonably–sized generalizable NN is not straightforward when
training data is limited to a narrow range in the shear–dominated region.
Local MAPE contours in the entire parameter space for DNNs with different
number of neurons are compared in Fig. 14. It is seen that the DNN with large
width have relatively smaller errors in both training and testing regions.

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

Figure 13:  MAPE of NNs with different architectures for Case–3 (a) training,
(b) testing

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(e)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(f)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(g)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(h)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(i)

Figure 14:  MAPE contours for Case–3, (a)-(c) 7L–5N, (d)-(f) 7L–7N, (g)-(i)
7L–15N

The main inferences from the three case studies can be summarized. We
systematically examined the suitability of DNNs for RANS turbulence closure
modeling. We studied the existence of a generalizable, moderate-sized NN given
sufficiently many training data in the entire parameter space. The results
show that NNs require more degrees of freedom (than the true proxy–physics
model) to accurately approximate the true polynomial model with four
coefficients even with entire data points of the parameter space. This implies
that DNNs are not efficient approximations of data manifolds with
non–linearity and bifurcation effects such as turbulence even in
interpolation. We also observed that the approximation capability of the NNs
significantly varies by activation function of neurons. Comparing the
performance of the NNs with different loss functions (MAPE and RMSE) showed
that networks with MAPE loss function have better performance for the
generated datasets in this study. For the interpolation case, it has been
shown that the reasonably–sized NNs trained with data manifold with no
bifurcation (SZL model) have the smallest errors. However, the coexistence of
non–linearity and bifurcation in the data manifold produced the largest level
of testing errors in the NN solutions even in these simple proxy–physics
models. As the non–linearity increases in the data manifolds with bifurcation,
the approximation capability of the NNs reduces in both interpolation and
extrapolation cases. Furthermore, in the extrapolation cases it is not
straightforward to find an optimal architecture.

##  VI  Conclusion

Generalizability of ML–assisted RANS turbulence model to unseen flows still
faces many challenges beck2021perspective  due to flow–dependent non–linearity
and bifurcations of the constitutive relations. Further, there is little
consensus and great deal of uncertainty regarding the choice of NN
hyperparameters and training techniques. Yet, these choices can significantly
affect the predictive capability and generalizability of ML turbulence models.
We seek to understand the optimal choice of hyperparameters, training process
elements (type of loss function) and necessary number of neurons of DNNs
required to allow a sufficiently accurate approximation at the RANS closure
modeling level. Standard fully–connected NNs are trained in a supervised
manner and their approximation capabilities are systematically investigated by
considering the effects of: _(i)_ intrinsic complexity of the solution
manifold; _(ii)_ sampling procedure (interpolation vs. extrapolation) and
_(iii)_ optimization procedure.

A key novelty of this work is the adoption of simplified proxy–physics
turbulence surrogates that incorporate some of the important features of real
homogeneous flows to generate the sufficient training data to assess
generalizability (interpolation vs. extrapolation) characteristics. An
important advantage of this approach is that training data for all flows in
the parameter space can be generated easily. In contrast, DNS or experiments
would be prohibitively expensive and may not even be feasible for all flows in
the parameter space. Successful generalizability of ML models in this
proxy–physics turbulence system is a necessary but not a sufficient condition
for ML–model generalizability in actual turbulent flows. Nevertheless, this
study provides valuable insight into the generalizability characteristics of
different network architectures and hyperparameters in turbulence–like
phenomena.

Three turbulence surrogates of different degrees of complexity are chosen:
_(i)_ a non–linear constitutive relation with no bifurcation in the parameter
regime _(ii)_ mildly non–linear constitutive relation with bifurcation in the
regime of interest; and _(iii)_ moderately non–linear constitutive model with
bifurcation. When the constitutive relation does not have bifurcation in the
entire parameter space, the accuracy of the ML model is quite reasonable.
However, when the surrogate data exhibits bifurcation, the combination of
non–linearity and bifurcation produced the largest level of testing errors
even in these simple proxy–physics models. Moreover, testing error increase
with the level of non–linearity in the proxy–physics surrogate with
bifurcation in the parameter space.

The conclusions of the study are two–fold. First, even for interpolation, the
NN–based models require very large networks (more degrees of freedom than true
proxy–physics model) to reduce errors to a reasonable level. Secondly, in a
practical calculation, it is difficult _a priori_ to guarantee that the model
will be queried only in the trained domain. Therefore, it is critical to
establish that the NN–model can handle certain level of extrapolation. To
assess the model ability for extrapolation/generalization, we train the model
in limited parameter space and test outside of this regime. It is shown that
even for this simple proxy–physics system, the NN–model performance is
inadequate (Cases 2, 3). Further, we identify and distinguish the challenges
to generalizability arising out of non–linearity and bifurcation. We believe
that these findings are important and yet not clearly understood in current
literature. Therefore, studies such as these are necessary for a balanced
assessment of NN–based RANS turbulence models as predictive tools for use in
unseen flows.

As mentioned before, true turbulence phenomena is much more complicated with
many more degrees of freedom (many more physical parameters) and multiple
bifurcations in the overall behavior. For example, mean flow
three–dimensionality, large–scale instabilities, streamline curvature,
stratification, compressibility and other effects encountered in typical
engineering flows will lead to significantly more complex constitutive
relations with non–equilibrium effects. It is very likely that sufficiently
large samples of accurate (direct numerical simulations or experimental) time
dependent data over all possible regimes of turbulence encompassing all these
effects will not be available in the near future. Even if such data were
available, it is unclear if an optimal set of hyperparameters can be found
using currently available methods. Until such a time when _(i)_ sufficiently
large volume of unsteady data is available over the entire parameter regime of
turbulence; and _(ii)_ analytical methods for determining optimal network
architecture are available, the results of the current study suggest that
ML–enhanced RANS models will be limited in its ability to perform predictive
computations of real engineering flows. It is important to note that when data
is indeed available over the entire parameter range, traditional closure
development methods may also lead to significantly improved models. Thus, the
relative advantages of ML methods over traditional methods must be reassessed
at that time.

##  APPENDIX

Table 5: Performance of models trained with different activation functions  Type |  MAPE–training |  MAPE–testing
---|---|---
ReLU |  0.93 |  8.53
Elu |  0.90 |  7.37
Leaky-ReLU |  0.5 |  11.60
Tanh |  0.75 |  12.03
Sigmoid |  0.85 |  4.39

For the Case–2 in which training data is partially available in the
strain–dominated region, a systematic grid search hyperparameter optimization
is performed. For this case, the data manifold is generated by ARSM with SSG
model. First, the performance of a 7L-7N NN with different activation
functions (act) is studied. During the training a $L_{2}$ norm regularization
with $\lambda=0.1$ is used. All other hyperparameters are fixed as shown in
Table 3. The training and testing MAPE errors for this experiment are shown in
Table 5. It is seen that between the five activation functions, ReLU, Elu, and
Sigmoid have reasonably small training and testing MAPE errors. By selecting
two ReLU and Sigmoid activation functions, the performance of the NNs with
different architectures are further investigated. MAPE of all the 16 network
architectures in the training and testing datasets is shown in Fig. 15.
Results show that for both activation functions, optimum network architecture
should be selected from the deep networks with small training and testing MAPE
errors. Additionally, it is seen that between two activation functions, larger
networks with Sigmoid activation function have smaller testing errors in this
case.

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(a)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(b)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(c)

\begin{picture}(0.0,0.0)\put(-138.0,0.0){}\end{picture}

(d)

Figure 15:  Training and testing MAPE for Case–2, (a) and (c) ReLU, (b) and
(d) Sigmoid

A 7L–7N NN with Sigmoid activation function are selected to study the effects
of other hyperparameters on training and performance of the models. The
results reported in Tables 6–9 justifies the selection of optimum
hyperparameters shown in Table 3 in this study.

Table 6: Performance of models trained with different learning rates, act=Sigmoid, bs=50, opt=Adam, $L_{2}(\lambda$=0.1) Learning rate |  MAPE–training |  MAPE–testing
---|---|---
$1.0\times 10^{-2}$ |  2.48 |  3.98
$1.0\times 10^{-3}$ |  0.85 |  4.39
$1.0\times 10^{-4}$ |  0.77 |  10.81
$1.0\times 10^{-5}$ |  0.74 |  10.98
Table 7: Performance of models trained with different batch sizes, act=Sigmoid, lr=0.001, opt=Adam, $L_{2}(\lambda$=0.1) Batch size |  MAPE–training |  MAPE–testing
---|---|---
25 |  0.88 |  4.41
50 |  0.85 |  4.39
100 |  0.91 |  9.63
1000 |  0.88 |  11.08
Table 8: Performance of models trained with different type of optimizers, act=Sigmoid, lr=0.001, bs=50, $L_{2}(\lambda$=0.1) Type |  MAPE–training |  MAPE–testing
---|---|---
Adam |  0.85 |  4.39
RMSProp |  1.92 |  5.43
Table 9: Performance of models trained with different regularization coefficients, act=Sigmoid, lr=0.001, bs=50 $L_{1}$–norm |  $L_{2}$–norm
---|---
Coefficient $\lambda$ |  Training |  Testing |  Training |  Testing
0 |  0.83 |  5.75 |  0.83 |  5.75
0.01 |  0.82 |  5.24 |  0.82 |  5.15
0.1 |  0.82 |  4.74 |  0.85 |  4.39
0.2 |  0.83 |  4.49 |  0.84 |  4.40

## Acknowledgement

The authors acknowledge support provided by Texas A&M High Performance
Research Computing center.

## References

## References

  * (1) G. Blaisdell and K. Shari, “Simulation and modeling of the elliptic streamline ow,” in Proceedings of the Summer Program, p. 1, Citeseer, 1996.
  * (2) S. S. Girimaji, “Pressure–strain correlation modelling of complex turbulent flows,” Journal of Fluid Mechanics, vol. 422, pp. 91–123, 2000.
  * (3) C. G. Speziale, “Analytical methods for the development of reynolds-stress closures in turbulence,” Annu. Rev. Fluid Mech, vol. 107, p. 57, 1991.
  * (4) A. A. Mishra and S. S. Girimaji, “Pressure–strain correlation modeling: towards achieving consistency with rapid distortion theory,” Flow, turbulence and combustion, vol. 85, no. 3, pp. 593–619, 2010.
  * (5) A. A. Mishra and S. S. Girimaji, “Toward approximating non-local dynamics in single-point pressure–strain correlation closures,” Journal of Fluid Mechanics, vol. 811, pp. 168–188, 2017.
  * (6) C. G. Speziale, R. Abid, and G. A. Blaisdell, “On the consistency of reynolds stress turbulence closures with hydrodynamic stability theory,” Physics of Fluids, vol. 8, no. 3, pp. 781–788, 1996.
  * (7) S. S. Girimaji, “Partially-averaged navier-stokes model for turbulence: A reynolds-averaged navier-stokes to direct numerical simulation bridging method,” Journal of applied mechanics, vol. 73, no. 3, pp. 413–421, 2006\.
  * (8) S. S. Girimaji, E. Jeong, and R. Srinivasan, “Partially averaged navier-stokes method for turbulence: Fixed point analysis and comparison with unsteady partially averaged navier-stokes,” Journal of Applied Mechanics, vol. 73, no. 3, p. 422, 2006.
  * (9) F. Sarghini, G. De Felice, and S. Santini, “Neural networks based subgrid scale modeling in large eddy simulations,” Computers & fluids, vol. 32, no. 1, pp. 97–108, 2003.
  * (10) M. Gamahara and Y. Hattori, “Searching for turbulence models by artificial neural network,” Physical Review Fluids, vol. 2, no. 5, p. 054604, 2017\.
  * (11) R. Maulik and O. San, “A neural network approach for the blind deconvolution of turbulent flows,” Journal of Fluid Mechanics, vol. 831, pp. 151–181, 2017.
  * (12) R. Maulik, O. San, A. Rasheed, and P. Vedula, “Data-driven deconvolution for large eddy simulations of kraichnan turbulence,” Physics of Fluids, vol. 30, no. 12, p. 125109, 2018.
  * (13) A. Beck, D. Flad, and C.-D. Munz, “Deep neural networks for data-driven les closure models,” Journal of Computational Physics, vol. 398, p. 108910, 2019.
  * (14) Z. Zhou, G. He, S. Wang, and G. Jin, “Subgrid-scale model for large-eddy simulation of isotropic turbulent flows using an artificial neural network,” Computers & Fluids, vol. 195, p. 104319, 2019.
  * (15) C. Xie, J. Wang, and E. Weinan, “Modeling subgrid-scale forces by spatial artificial neural networks in large eddy simulation of turbulence,” Physical Review Fluids, vol. 5, no. 5, p. 054606, 2020.
  * (16) M. Schoepplein, J. Weatheritt, R. Sandberg, M. Talei, and M. Klein, “Application of an evolutionary algorithm to les modelling of turbulent transport in premixed flames,” Journal of Computational Physics, vol. 374, pp. 1166–1179, 2018.
  * (17) M. Reissmann, J. Hasslberger, R. D. Sandberg, and M. Klein, “Application of gene expression programming to a-posteriori les modeling of a taylor green vortex,” Journal of Computational Physics, vol. 424, p. 109859, 2021.
  * (18) R. Maulik, O. San, A. Rasheed, and P. Vedula, “Subgrid modelling for two-dimensional turbulence using neural networks,” Journal of Fluid Mechanics, vol. 858, pp. 122–144, 2019.
  * (19) Z. Wang, K. Luo, D. Li, J. Tan, and J. Fan, “Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation,” Physics of Fluids, vol. 30, no. 12, p. 125101, 2018.
  * (20) R. Maulik, O. San, J. D. Jacob, and C. Crick, “Sub-grid scale model classification and blending through deep learning,” Journal of Fluid Mechanics, vol. 870, pp. 784–812, 2019.
  * (21) C. Xie, J. Wang, H. Li, M. Wan, and S. Chen, “Artificial neural network mixed model for large eddy simulation of compressible isotropic turbulence,” Physics of Fluids, vol. 31, no. 8, p. 085112, 2019.
  * (22) C. Xie, Z. Yuan, and J. Wang, “Artificial neural network-based nonlinear algebraic models for large eddy simulation of turbulence,” Physics of Fluids, vol. 32, no. 11, p. 115101, 2020.
  * (23) C. Xie, K. Li, C. Ma, and J. Wang, “Modeling subgrid-scale force and divergence of heat flux of compressible isotropic turbulence by artificial neural network,” Physical Review Fluids, vol. 4, no. 10, p. 104605, 2019\.
  * (24) J. Park and H. Choi, “Toward neural-network-based large eddy simulation: application to turbulent channel flow,” Journal of Fluid Mechanics, vol. 914, 2021.
  * (25) Y. Wang, Z. Yuan, C. Xie, and J. Wang, “Artificial neural network-based spatial gradient models for large-eddy simulation of turbulence,” AIP Advances, vol. 11, no. 5, p. 055216, 2021.
  * (26) A. Subel, A. Chattopadhyay, Y. Guan, and P. Hassanzadeh, “Data-driven subgrid-scale modeling of forced burgers turbulence using deep learning with generalization to higher reynolds numbers via transfer learning,” Physics of Fluids, vol. 33, no. 3, p. 031702, 2021.
  * (27) S. Pawar, O. San, A. Rasheed, and P. Vedula, “A priori analysis on deep learning of subgrid-scale parameterizations for kraichnan turbulence,” Theoretical and Computational Fluid Dynamics, vol. 34, no. 4, pp. 429–455, 2020\.
  * (28) J. Ling, A. Kurzawski, and J. Templeton, “Reynolds averaged turbulence modelling using deep neural networks with embedded invariance,” Journal of Fluid Mechanics, vol. 807, pp. 155–166, 2016.
  * (29) Z. Zhang, X.-d. Song, S.-r. Ye, Y.-w. Wang, C.-g. Huang, Y.-r. An, and Y.-s. Chen, “Application of deep learning method to reynolds stress models of channel flow based on reduced-order modeling of dns data,” Journal of Hydrodynamics, vol. 31, no. 1, pp. 58–65, 2019.
  * (30) C. Sotgiu, B. Weigand, K. Semmler, and P. Wellinger, “Towards a general data-driven explicit algebraic reynolds stress prediction framework,” International Journal of Heat and Fluid Flow, vol. 79, p. 108454, 2019.
  * (31) N. Geneva and N. Zabaras, “Quantifying model form uncertainty in reynolds-averaged turbulence models with bayesian deep neural networks,” Journal of Computational Physics, vol. 383, pp. 125–147, 2019.
  * (32) R. Fang, D. Sondak, P. Protopapas, and S. Succi, “Neural network models for the anisotropic reynolds stress tensor in turbulent channel flow,” Journal of Turbulence, vol. 21, no. 9-10, pp. 525–543, 2020.
  * (33) C. Jiang, J. Mi, S. Laima, and H. Li, “A novel algebraic stress model with machine-learning-assisted parameterization,” Energies, vol. 13, no. 1, p. 258, 2020.
  * (34) Y. Yin, P. Yang, Y. Zhang, H. Chen, and S. Fu, “Feature selection and processing of turbulence modeling based on an artificial neural network,” Physics of Fluids, vol. 32, no. 10, p. 105117, 2020.
  * (35) S. Taghizadeh, F. D. Witherden, and S. S. Girimaji, “Turbulence closure modeling with data-driven techniques: physical compatibility and consistency considerations,” New Journal of Physics, vol. 22, no. 9, p. 093023, 2020\.
  * (36) C. Jiang, R. Vinuesa, R. Chen, J. Mi, S. Laima, and H. Li, “An interpretable framework of data-driven turbulence modeling using deep neural networks,” Physics of Fluids, vol. 33, no. 5, p. 055133, 2021.
  * (37) X.-H. Zhou, J. Han, and H. Xiao, “Learning nonlocal constitutive models with neural networks,” Computer Methods in Applied Mechanics and Engineering, vol. 384, p. 113927, 2021.
  * (38) M. Yang and Z. Xiao, “Improving the k–$\omega$–$\gamma$–ar transition model by the field inversion and machine learning framework,” Physics of Fluids, vol. 32, no. 6, p. 064101, 2020.
  * (39) P. S. Volpiani, M. Meyer, L. Franceschini, J. Dandois, F. Renac, E. Martin, O. Marquet, and D. Sipp, “Machine learning-augmented turbulence modeling for rans simulations of massively separated flows,” Physical Review Fluids, vol. 6, no. 6, p. 064607, 2021.
  * (40) R. Maulik, H. Sharma, S. Patel, B. Lusch, and E. Jennings, “A turbulent eddy-viscosity surrogate modeling framework for reynolds-averaged navier-stokes simulations,” Computers & Fluids, vol. 227, p. 104777, 2021\.
  * (41) L. Zhu, W. Zhang, J. Kou, and Y. Liu, “Machine learning methods for turbulence modeling in subsonic flows around airfoils,” Physics of Fluids, vol. 31, no. 1, p. 015105, 2019.
  * (42) E. L. Peters, R. Balin, K. E. Jansen, A. Doostan, and J. A. Evans, “S-frame discrepancy correction models for data-informed reynolds stress closure,” Journal of Computational Physics, p. 110717, 2021.
  * (43) L. Breiman, “Random forests,” Machine learning, vol. 45, no. 1, pp. 5–32, 2001.
  * (44) J.-L. Wu, J.-X. Wang, H. Xiao, and J. Ling, “A priori assessment of prediction confidence for data-driven turbulence modeling,” Flow, Turbulence and Combustion, vol. 99, no. 1, pp. 25–46, 2017.
  * (45) M. L. Kaandorp and R. P. Dwight, “Data-driven modelling of the reynolds stress tensor using random forests with invariance,” Computers & Fluids, vol. 202, p. 104497, 2020.
  * (46) J.-L. Wu, H. Xiao, and E. Paterson, “Physics-informed machine learning approach for augmenting turbulence models: A comprehensive framework,” Physical Review Fluids, vol. 3, no. 7, p. 074602, 2018.
  * (47) J. Weatheritt and R. Sandberg, “The development of algebraic stress models using a novel evolutionary algorithm,” International Journal of Heat and Fluid Flow, vol. 68, pp. 298–318, 2017.
  * (48) J. Weatheritt, R. Pichler, R. D. Sandberg, G. Laskowski, and V. Michelassi, “Machine learning for turbulence model development using a high-fidelity hpt cascade simulation,” in ASME Turbo Expo 2017: Turbomachinery Technical Conference and Exposition, American Society of Mechanical Engineers Digital Collection, 2017.
  * (49) H. D. Akolekar, J. Weatheritt, N. Hutchins, R. D. Sandberg, G. Laskowski, and V. Michelassi, “Development and use of machine-learnt algebraic reynolds stress models for enhanced prediction of wake mixing in lpts,” in ASME Turbo Expo 2018: Turbomachinery Technical Conference and Exposition, American Society of Mechanical Engineers Digital Collection, 2018.
  * (50) Y. Zhao, H. D. Akolekar, J. Weatheritt, V. Michelassi, and R. D. Sandberg, “Rans turbulence model development using cfd-driven machine learning,” Journal of Computational Physics, vol. 411, p. 109413, 2020.
  * (51) C. Lav, R. D. Sandberg, and J. Philip, “A framework to develop data-driven turbulence models for flows with organised unsteadiness,” Journal of Computational Physics, vol. 383, pp. 148–165, 2019.
  * (52) F. Waschkowski, Y. Zhao, R. Sandberg, and J. Klewicki, “Multi-objective cfd-driven development of coupled turbulence closure models,” arXiv preprint arXiv:2105.06225, 2021.
  * (53) M. Schmelzer, R. P. Dwight, and P. Cinnella, “Discovery of algebraic reynolds-stress models using sparse symbolic regression,” Flow, Turbulence and Combustion, vol. 104, no. 2, pp. 579–603, 2020.
  * (54) J. P. Huijing, R. P. Dwight, and M. Schmelzer, “Data-driven rans closures for three-dimensional flows around bluff bodies,” Computers & Fluids, vol. 225, p. 104997, 2021.
  * (55) S. Beetham and J. Capecelatro, “Formulating turbulence closures using sparse regression with embedded form invariance,” Physical Review Fluids, vol. 5, no. 8, p. 084611, 2020.
  * (56) Y. Zhang, R. P. Dwight, M. Schmelzer, J. F. Gómez, Z.-h. Han, and S. Hickel, “Customized data-driven rans closures for bi-fidelity les–rans optimization,” Journal of Computational Physics, vol. 432, p. 110153, 2021\.
  * (57) J. N. Kutz, “Deep learning in fluid dynamics,” Journal of Fluid Mechanics, vol. 814, pp. 1–4, 2017.
  * (58) K. Duraisamy, G. Iaccarino, and H. Xiao, “Turbulence modeling in the age of data,” Annual Review of Fluid Mechanics, vol. 51, pp. 357–377, 2019.
  * (59) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine learning for fluid mechanics,” Annual Review of Fluid Mechanics, vol. 52, pp. 477–508, 2020\.
  * (60) A. Beck and M. Kurz, “A perspective on machine learning methods in turbulence modeling,” GAMM-Mitteilungen, vol. 44, no. 1, p. e202100002, 2021.
  * (61) K. Duraisamy, “Perspectives on machine learning-augmented reynolds-averaged and large eddy simulation models of turbulence,” Physical Review Fluids, vol. 6, no. 5, p. 050504, 2021.
  * (62) K. Anand, Z. Wang, M. Loog, and J. van Gemert, “Black magic in deep learning: How human skill impacts network training,” arXiv preprint arXiv:2008.05981, 2020.
  * (63) G. Cybenko, “Approximation by superpositions of a sigmoidal function,” Mathematics of control, signals and systems, vol. 2, no. 4, pp. 303–314, 1989\.
  * (64) K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural networks, vol. 2, no. 5, pp. 359–366, 1989.
  * (65) A. R. Barron, “Universal approximation bounds for superpositions of a sigmoidal function,” IEEE Transactions on Information theory, vol. 39, no. 3, pp. 930–945, 1993.
  * (66) M. Leshno, V. Y. Lin, A. Pinkus, and S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function,” Neural networks, vol. 6, no. 6, pp. 861–867, 1993.
  * (67) V. Maiorov and A. Pinkus, “Lower bounds for approximation by mlp neural networks,” Neurocomputing, vol. 25, no. 1-3, pp. 81–91, 1999.
  * (68) H. N. Mhaskar, “Approximation properties of a multilayered feedforward artificial neural network,” Advances in Computational Mathematics, vol. 1, no. 1, pp. 61–80, 1993.
  * (69) H. N. Mhaskar, “Neural networks for optimal approximation of smooth and analytic functions,” Neural computation, vol. 8, no. 1, pp. 164–177, 1996\.
  * (70) D. Yarotsky, “Error bounds for approximations with deep relu networks,” Neural Networks, vol. 94, pp. 103–114, 2017.
  * (71) P. Petersen and F. Voigtlaender, “Optimal approximation of piecewise smooth functions using deep relu neural networks,” Neural Networks, vol. 108, pp. 296–330, 2018.
  * (72) M. Raissi, P. Perdikaris, and G. E. Karniadakis, “Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations,” Journal of Computational Physics, vol. 378, pp. 686–707, 2019.
  * (73) X. Jin, S. Cai, H. Li, and G. E. Karniadakis, “Nsfnets (navier-stokes flow nets): Physics-informed neural networks for the incompressible navier-stokes equations,” Journal of Computational Physics, vol. 426, p. 109951, 2021\.
  * (74) J. Berner, P. Grohs, and A. Jentzen, “Analysis of the generalization error: Empirical risk minimization over deep artificial neural networks overcomes the curse of dimensionality in the numerical approximation of black–scholes partial differential equations,” SIAM Journal on Mathematics of Data Science, vol. 2, no. 3, pp. 631–657, 2020.
  * (75) E. Weinan, J. Han, and A. Jentzen, “Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations,” Communications in Mathematics and Statistics, vol. 5, no. 4, pp. 349–380, 2017.
  * (76) J. Han, A. Jentzen, and E. Weinan, “Solving high-dimensional partial differential equations using deep learning,” Proceedings of the National Academy of Sciences, vol. 115, no. 34, pp. 8505–8510, 2018.
  * (77) J. Sirignano and K. Spiliopoulos, “Dgm: A deep learning algorithm for solving partial differential equations,” Journal of computational physics, vol. 375, pp. 1339–1364, 2018.
  * (78) N. Dal Santo, S. Deparis, and L. Pegolotti, “Data driven approximation of parametrized pdes by reduced basis and neural networks,” Journal of Computational Physics, vol. 416, p. 109550, 2020.
  * (79) J. S. Hesthaven and S. Ubbiali, “Non-intrusive reduced order modeling of nonlinear problems using neural networks,” Journal of Computational Physics, vol. 363, pp. 55–78, 2018.
  * (80) Y. Khoo, J. Lu, and L. Ying, “Solving parametric pde problems with artificial neural networks,” European Journal of Applied Mathematics, vol. 32, no. 3, pp. 421–435, 2021.
  * (81) K. Lee and K. T. Carlberg, “Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders,” Journal of Computational Physics, vol. 404, p. 108973, 2020.
  * (82) G. Kutyniok, P. Petersen, M. Raslan, and R. Schneider, “A theoretical analysis of deep neural networks and parametric pdes,” Constructive Approximation, pp. 1–53, 2021.
  * (83) M. Geist, P. Petersen, M. Raslan, R. Schneider, and G. Kutyniok, “Numerical solution of the parametric diffusion equation by deep neural networks,” Journal of Scientific Computing, vol. 88, no. 1, pp. 1–37, 2021.
  * (84) S. S. Girimaji, “Fully explicit and self-consistent algebraic reynolds stress model,” Theoretical and Computational Fluid Dynamics, vol. 8, no. 6, pp. 387–402, 1996.
  * (85) T. Gatski, “On explicit algebraic stress models for complex turbulent flows,” J. Fluid Mech, vol. 254, pp. 59–78, 1993.
  * (86) B. E. Launder, G. J. Reece, and W. Rodi, “Progress in the development of a reynolds-stress turbulence closure,” Journal of fluid mechanics, vol. 68, no. 3, pp. 537–566, 1975.
  * (87) C. SPEZIALE, S. SARKAR, and T. GATSKI, “Modelling the pressure-strain correlation of turbulence: an invariant dynamical stystems approach,” Journal of Fluid Mechanics, vol. 227, pp. 245–272, 1991.
  * (88) S. S. Girimaji, “Lower-dimensional manifold (algebraic) representation of reynolds stress closure equations,” Theoretical and Computational Fluid Dynamics, vol. 14, no. 4, pp. 259–281, 2001.
  * (89) C. A. Gomez and S. S. Girimaji, “Explicit algebraic reynolds stress model (earsm) for compressible shear flows,” Theoretical and Computational Fluid Dynamics, vol. 28, no. 2, pp. 171–196, 2014.
  * (90) A. A. Mishra and S. S. Girimaji, “Intercomponent energy transfer in incompressible homogeneous turbulence: multi-point physics and amenability to one-point closures,” Journal of Fluid Mechanics, vol. 731, pp. 639–681, 2013.
  * (91) T.-H. Shih, J. Zhu, and J. L. Lumley, A realizable Reynolds stress algebraic equation model, vol. 105993.  National Aeronautics and Space Administration, 1993.
  * (92) T. Craft, B. Launder, and K. Suga, “Development and application of a cubic eddy-viscosity model of turbulence,” International Journal of Heat and Fluid Flow, vol. 17, no. 2, pp. 108–115, 1996.
  * (93) I. Goodfellow, Y. Bengio, and A. Courville, Deep learning, vol. 1.  MIT press Cambridge, 2016.
  * (94) S. Cai, Z. Wang, L. Lu, T. A. Zaki, and G. E. Karniadakis, “Deepm&mnet: Inferring the electroconvection multiphysics fields based on operator approximation by neural networks,” Journal of Computational Physics, vol. 436, p. 110296, 2021.
  * (95) G. James, D. Witten, T. Hastie, and R. Tibshirani, An introduction to statistical learning, vol. 112.  Springer, 2013.
  * (96) D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” arXiv preprint arXiv:1412.6980, 2014.
  * (97) X. Glorot and Y. Bengio, “Understanding the difficulty of training deep feedforward neural networks,” in Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249–256, JMLR Workshop and Conference Proceedings, 2010.