Connecting implicit and explicit large eddy simulations of two-dimensional turbulence through machine learning

TL;DR

This paper introduces a machine learning framework that dynamically chooses between implicit and explicit dissipation schemes in large eddy simulations of 2D turbulence, improving statistical accuracy without predefined turbulence models.

Contribution

It presents a novel data-driven approach to adaptively select numerical dissipation schemes in LES based on DNS data, bridging implicit and explicit turbulence modeling.

Findings

Enhanced statistical fidelity in LES results.

Effective dynamic scheme selection based on learned criteria.

Improved energy spectra and vorticity structures.

Abstract

In this article, we utilize machine learning to dynamically determine if a point on the computational grid requires implicit numerical dissipation for large eddy simulation (LES). The decision making process is learnt through \emph{a priori} training on quantities derived from direct numerical simulation (DNS) data. In particular, we compute eddy-viscosities obtained through the coarse graining of DNS quantities and utilize their distribution to categorize areas that require dissipation. If our learning determines that closure is necessary, an upwinded scheme is utilized for computing the non-linear Jacobian. In contrast, if it is determined that closure is unnecessary, a symmetric and second-order accurate energy and enstrophy preserving Arakawa scheme is utilized instead. This results in a closure framework that precludes the specification of any model-form for the small scale contributions of turbulence but deploys an appropriate numerical dissipation from explicit closure driven hypotheses. This methodology is deployed for the Kraichnan turbulence test-case and assessed through various statistical quantities such as angle-averaged kinetic energy spectra and vorticity structure functions. Our framework thus establishes a direct link between the use of explicit LES ideologies for closure and numerical scheme-based modeling of turbulence leading to improved statistical fidelity of \emph{a posteriori} simulations.

Figures

Equations

$$\frac{\partial \omega }{\partial t}+J(\omega ,\psi )=\frac{1}{Re}{\nabla }^2\omega ,$$

$${\nabla }^2\psi =-\omega .$$

$$J(\omega ,\psi )=\frac{\partial \psi }{\partial y}\frac{\partial \omega }{\partial x}-\frac{\partial \psi }{\partial x}\frac{\partial \omega }{\partial y}.$$

$$u$$

$$\frac{\partial \bar{\omega }}{\partial t}+J(\bar{\omega },\bar{\psi })=\frac{1}{Re}{\nabla }^2\bar{\omega },$$

$$J^A(\bar{\omega },\bar{\psi })=\frac{1}{3}\left(J_1(\bar{\omega },\bar{\psi })+J_2(\bar{\omega },\bar{\psi })+J_3(\bar{\omega },\bar{\psi })\right)$$

$$\begin{array}{rl}{J}_1(\bar{\omega },\bar{\psi })& =\frac{1}{4\Delta x\Delta y}[({\bar{\omega }}_{i+1,j}-{\bar{\omega }}_{i-1,j})({\bar{\psi }}_{i,j+1}-{\bar{\psi }}_{i,j-1})\\ & -({\bar{\omega }}_{i,j+1}-{\bar{\omega }}_{i,j-1})({\bar{\psi }}_{i+1,j}-{\bar{\psi }}_{i-1,j})],\end{array}$$

$$\begin{array}{rl}& J_2(\bar{\omega },\bar{\psi })=\frac{1}{4\Delta x\Delta y}[{\bar{\omega }}_{i+1,j}({\bar{\psi }}_{i+1,j+1}-{\bar{\psi }}_{i+1,j-1})\\ & -{\bar{\omega }}_{i-1,j}({\bar{\psi }}_{i-1,j+1}-{\bar{\psi }}_{i-1,j-1})-{\bar{\omega }}_{i,j+1}({\bar{\psi }}_{i+1,j+1}-{\bar{\psi }}_{i-1,j+1})\\ & +{\bar{\omega }}_{i,j-1}({\bar{\psi }}_{i+1,j-1}-{\bar{\psi }}_{i-1,j-1})],\end{array}$$

$$\begin{array}{rl}& J_3(\bar{\omega },\bar{\psi })=\frac{1}{4\Delta x\Delta y}[{\bar{\omega }}_{i+1,j+1}({\bar{\psi }}_{i,j+1}-{\bar{\psi }}_{i+1,j})\\ & -{\bar{\omega }}_{i-1,j-1}({\bar{\psi }}_{i-1,j}-{\bar{\psi }}_{i,j-1})-{\bar{\omega }}_{i-1,j+1}({\bar{\psi }}_{i,j+1}-{\bar{\psi }}_{i-1,j})\\ & +{\bar{\omega }}_{i+1,j-1}({\bar{\psi }}_{i+1,j}-{\bar{\psi }}_{i,j-1})].\end{array}$$

$$\begin{array}{rl}{J}^I(\bar{\omega },\bar{\psi })=& {\bar{u}}_{i,j}\frac{{\bar{\omega }}_{i+1,j}-{\bar{\omega }}_{i-1,j}}{2\Delta x}+\frac{1}{2}({\bar{u}}^{+}{\bar{\omega }}_x^{-}+u^{-}{\bar{\omega }}_x^{+})\\ & +{\bar{v}}_{i,j}\frac{{\bar{\omega }}_{i,j+1}-{\bar{\omega }}_{i,j-1}}{2\Delta y}+\frac{1}{2}({\bar{v}}^{+}{\bar{\omega }}_y^{-}+{\bar{v}}^{-}{\bar{\omega }}_y^{+}),\end{array}$$

$${\bar{u}}^{-}$$

$${\bar{v}}^{-}$$

$$\begin{array}{rl}{\bar{\omega }}_x^{-}& =\frac{{\bar{\omega }}_{i-2,j}-3{\bar{\omega }}_{i-1,j}+3{\bar{\omega }}_{i,j}-{\bar{\omega }}_{i+1,j}}{3\Delta x},\\ {\bar{\omega }}_x^{+}& =\frac{{\bar{\omega }}_{i-1,j}-3{\bar{\omega }}_{i,j}+3{\bar{\omega }}_{i+1,j}-{\bar{\omega }}_{i+2,j}}{3\Delta x},\\ {\bar{\omega }}_y^{-}& =\frac{{\bar{\omega }}_{i,j-2}-3{\bar{\omega }}_{i,j-1}+3{\bar{\omega }}_{i,j}-{\bar{\omega }}_{i,j+1}}{3\Delta y},\\ {\bar{\omega }}_y^{+}& =\frac{{\bar{\omega }}_{i,j-1}-3{\bar{\omega }}_{i,j}+3{\bar{\omega }}_{i,j+1}-{\bar{\omega }}_{i,j+2}}{3\Delta y}.\end{array}$$

$$\begin{array}{rl}{\bar{u}}_{i,j}& =\frac{{\bar{\psi }}_{i,j+1}-{\bar{\psi }}_{i,j-1}}{2\Delta y}\\ {\bar{v}}_{i,j}& =-\frac{{\bar{\psi }}_{i+1,j}-{\bar{\psi }}_{i-1,j}}{2\Delta x},\end{array}$$

$$\Pi =J(\bar{\omega },\bar{\psi })-\overline{J(\omega ,\psi )}$$

$${\nu }_e^a=\frac{\Pi }{{\nabla }^2\bar{\omega }}$$

$$\mathbb{M}:\{{\mathcal{P}}_1,{\mathcal{P}}_2,\dots ,{\mathcal{P}}_P\}\in {\mathbb{R}}^P\to \{{\mathcal{Q}}_1,{\mathcal{Q}}_2,\dots ,{\mathcal{Q}}_Q\}\in {\mathbb{R}}^Q.$$

$$\mathbb{M}:\{\text{p}\}\in {\mathbb{R}}^{18}\to \{P(\text{q}|\text{p})\}\in {\mathbb{R}}^3,$$

$$\begin{array}{c}{\text{p}}_{i,j}=\{{\bar{\omega }}_{i,j},{\bar{\omega }}_{i,j+1},{\bar{\omega }}_{i,j-1},\dots ,{\bar{\omega }}_{i-1,j-1},\\ {\bar{\psi }}_{i,j},{\bar{\psi }}_{i,j+1},{\bar{\psi }}_{i,j-1},\dots ,{\bar{\psi }}_{i-1,j-1}\}\end{array}$$

$$P(\text{q}|\text{p}{)}_{i,j}=\{P(J^k(\bar{\omega },\bar{\psi }{)}_{i,j}|{\text{p}}_{i,j})\},$$

$$E(\text{w})=-\sum _{n=1}^N\sum _{k=1}^K\{t_{nk}\log (y_{nk})+(1-t_{nk})\log (1-y_{nk})\},$$

$$S_{\omega }^x=⟨|\omega (x+r,y)-\omega (x,y){|}^2⟩$$

$$S_{\omega }^y=⟨|\omega (x,y+r)-\omega (x,y){|}^2⟩,$$

Full text

Turbulence modeling, Machine learning

I Introduction

Over the past decade, advances in data collection and increasing access to computational resources have led to a revolution in the use of data-driven techniques for the solution of complex inverse problems. One such problem is that of turbulence, the multiscale nature of which causes extreme computational demands for most practical systems. As a result, turbulence requires the use multiple modeling approximations for the higher wavenumbers which remain unsupported by computational degrees of freedom. One such modeling approach is that of large eddy simulation (LES) Sagaut (2006), which attempts to simulate the evolution of the smaller wavenumbers while the unresolved frequencies are modeled by an algebraic or differential equation. As such, the basic premise of LES is extendable to many partial differential equation systems with quadratic non-linearities. The procedure of modeling these smaller scales is often denoted closure due to insufficient knowledge about higher-order wavenumber interactions with the coarse-grained system Berselli, Iliescu, and Layton (2006) and remains vital for the accurate computation of many applications Hickel, Egerer, and Larsson (2014); Yu, Xiao, and Li (2016); Zhou, Wang, and Jin (2018). From an LES point of view, the closure problem may be considered to be dominated by commutative errors in the calculation of the non-linear term as well as the defects associated with commutative errors stemming from the dynamic term. In this study, we focus on the former.

There are two main schools of thought when it comes to the LES of the Navier-Stokes equations. The first of these promotes the use of explicit closures. Explicit LES argues for the utilization of closures in the form of sub-grid models specified as algebraic or differential equations for the unresolved scales. These are built on intuitive reasoning of how the losses of coarse graining the Navier-Stokes equations may be incorporated into an LES deployment. Some of the most notable sub-grid closure strategies are those given by the eddy-viscosity hypothesis. Within the context of the Navier-Stokes equations, it is generally accepted that the finer scales are dissipative at the Kolmogorov length scales Kolmogorov (1941) and therefore, most turbulence models seek to specify a sub-grid dissipation Frisch (1995). Most sub-grid models can be traced back to the seminal work of Smagorinsky Smagorinsky (1963), where a model was proposed based on the concepts of an effective eddy-viscosity determined by an a priori specified mixing length and a $k^{-5/3}$ scaling recovery for the kinetic energy content in the wavenumber domain. Similar hypotheses have also been used for two-dimensional turbulence Leith (1968) (often utilized as a test-bed for geophysical scenarios, for instance see works by Pearson et al.Pearson and Fox-Kemper (2018); Pearson et al. (2017)), for approximating the $k^{-3}$ cascade in two-dimensional turbulence and generally have their roots in dimensional analysis related to the cascade of enstrophy. These models may also be classified as functional due to the phenomenological nature of their deployment and represent the bulk of explicit LES turbulence models used in practical deployments. Explicit LES closures may also be specified through the specification of a low-pass spatial filter to account for the unresolved scales Bardina, Ferziger, and Reynolds (1980); Stolz and Adams (1999); Layton and Lewandowski (2003); Mathew et al. (2003) where phenomenology is bypassed but ansatz are provided for the bulk dissipative nature of the smaller scales through the control of a characteristic filter-width. In either scenario, (i.e., whether structural or functional), the choice of the phenomenology (or dissipation control parameter) plays a key role in the successful calculation of accurate a posteriori statistics. In contrast, the implicit LES (or ILES) approach utilizes numerical dissipation to model the unresolved scales in a turbulent flow Grinstein, Margolin, and Rider (2007); El Rafei, Könözsy, and Rana (2017); Margolin (2018). In essence, the predominantly dissipative effects of the smallest scales are replicated through an artificial numerical dissipation via a biased discretization used in the calculation of the non-linear advective term Thornber, Mosedale, and Drikakis (2007); DeBonis (2013). The ILES approach is popular due to reduced algorithmic complexity and represents a union of turbulence modeling and shock capturing mechanisms but is often criticized due to the difficulties involved in quantifying the correct amount of dissipation in a turbulent flow evolution. This results in ILES methods often proving robust and stable but overly dissipative. In this work, we propose a machine learning algorithm to enable selective dissipation within an ILES deployment through the use of explicit LES concepts during the training of the learning framework.

The past few years have seen a rapid increase in the use of data-driven techniques for the spatio-temporal modeling of dynamical systems Schmidt and Lipson (2009); Bright, Lin, and Kutz (2013); Xiao et al. (2015); Ma, Lu, and Tryggvason (2015); Gautier et al. (2015); Brunton, Proctor, and Kutz (2016); Schaeffer (2017); Raissi, Perdikaris, and Karniadakis (2017); Mohan and Gaitonde (2018); Raissi and Karniadakis (2018); Rudy, Kutz, and Brunton (2018); San and Maulik (2018); Wan et al. (2018); Kim et al. (2018); Muravleva, Oseledets, and Koroteev (2018); Jin et al. (2018). When it comes to turbulence, some widely used strategies for inference include symbolic regression Weatheritt and Sandberg (2016, 2017a, 2017b), where functional model-forms for Reynolds-averaged Navier-Stokes (RANS) deployments were generated through evolutionary optimization against high-fidelity data. Other techniques incorporating Bayesian ideologies have also been used, for instance by Xiao et al.Xiao et al. (2016) where an iterative ensemble Kalman method was used to assimilate prior data for quantifying model form uncertainty in RANS models. In Wang et al.Wang, Wu, and Xiao (2017); Wang et al. (2017) and Wu et al.Wu, Xiao, and Paterson (2018), random-forest regressors were utilized for RANS turbulence-modeling given direct numerical simulation (DNS) data. In Singh and Duraisamy Singh and Duraisamy (2016) and Singh et al.Singh, Medida, and Duraisamy (2017), an ANN was utilized to predict a non-dimensional correction factor in the Spalart-Allmaras turbulence model through a field-inversion process using experimental data. Bypassing functional formulations of a turbulence model (a focus of this study) was also studied from the RANS point of view by Tracey et al. Tracey, Duraisamy, and Alonso (2015). Ling and Templeton Ling and Templeton (2015) utilized support vector machines, decision trees and random forest regressors for identifying regions of high RANS uncertainty. A deep-learning framework where Reynolds-stresses would be predicted in an invariant subspace was developed by Ling et al. Ling, Kurzawski, and Templeton (2016). Machine learning of invariance properties has also been discussed in the context of turbulence modeling Ling, Jones, and Templeton (2016). The reader is directed to a recent review by Duraisamy et al.Duraisamy, Iaccarino, and Xiao (2018), for an excellent review of turbulence modeling using data-driven ideas.

As shown above, the use of data-driven ideologies and in particular artificial neural networks (ANNs) has generated significant interest in the turbulence modeling community for addressing long-standing challenges (also see Sotgiu, Weigand, and Semmler (2018); Zhang et al. (2018); Zhu et al. (2019); Zhang et al. (2019); Raissi et al. (2019) for recent progress). One motivation for the popularity of ANNs is that a multilayered ANN may be optimally trained to universally approximate any non-linear function Hornik, Stinchcombe, and White (1989). In addition, the deployment of ANNs is amenable to integration within existing computational frameworks. Greater accessibility to data and ever-improving computing capabilities has also motivated the development of advanced ANN architectures for large-scale learning of complicated physical phenomena such as turbulence. Within the context of LES (and associated with the scope of this paper) there are several investigations into sub-grid modeling using data-driven techniques. In one of the first studies of the feasibility of using learning from DNS based high-fidelity data, Sarghini et al.Sarghini, De Felice, and Santini (2003) utilized ANNs for estimating Smagorinsky model-form coefficients within a mixed sub-grid model for a turbulent channel flow. ANNs were also used for wall-modeling by Milano and Koumotsakos Milano and Koumoutsakos (2002) where it was used to reconstruct the near wall field and compared to standard proper-orthogonal-decomposition techniques. An alternative to ANNs for sub-grid predictions was proposed by King et al.King, Hamlington, and Dahm (2016) where a priori optimization was utilized to minimize the $L^{2}$-error between true and modeled sub-grid quantities in a least-squares sense using a parameter-free Volterra series. Maulik and San Maulik and San (2017a) utilized an extreme-learning-machine (a variant of a single-layered ANN) to obtain maps between low-pass spatially filtered and deconvolved variables in an a priori sense. This had implications for the use of ANNs for turbulence modeling without model-form specification. A more in-depth investigation was recently undertaken by Fukami et al.Fukami, Fukagata, and Taira (2018) where convolutional ANNs were utilized for reconstructing from downsampled snapshots of turbulence. Maulik et al. Maulik et al. (2018a) also deployed a data-driven convolutional and deconvolutional operation to obtain closure terms for two-dimensional turbulence. Gamahara and Hattori Gamahara and Hattori (2017), utilized ANNs for identifying correlations with grid-resolved quantities for an indirect method of model-form identification in turbulent channel flow. The study by Vollant et al. Vollant, Balarac, and Corre (2017) utilized ANNs in conjuction with optimal estimator theory to obtain functional forms for sub-grid stresses. In Beck et al.Beck, Flad, and Munz (2018), a variety of neural network architectures such as convolutional and recurrent neural networks are studied for predicting closure terms for decaying homogeneous isotropic turbulence. A least-squares based truncation is specified for stable deployments of their model-free closures. Model-free turbulence closures are also specified by Maulik et al.Maulik et al. (2018a, 2019) and Wang et al.Wang et al. (2018), where sub-grid scale stresses are learned directly from DNS data and deployed in a posteriori assessments. King et al.King et al. (2018) studied generative-adversarial networks and the LAT-NET Hennigh (2017) for a priori recovery of statistics such as the intermittency of turbulent fluctuations and spectral scaling. A detailed discussion of the potential benefits and challenges of deep learning for turbulence (and fluid dynamics in general) may be found in the article by Kutz Kutz (2017).

While a large majority of the LES-based frameworks presented above utilize a least-squares error minimization technique for constructing maps to sub-grid stresses directly for theoretically optimal LES Langford and Moser (1999); Moser et al. (2009); LaBryer, Attar, and Vedula (2015), this work is novel in that it utilizes sub-grid statistics (pre-computed from DNS data) to train a classifier. This classifier determines whether a location requires dissipation or not through a priori experience in the learning phase. Once classified, the non-linear term at this particular point is evaluated using one of two schemes. If it is determined that the point requires no sub-grid closure, a symmetric and second-order accurate, energy and enstrophy conserving Arakawa-scheme Arakawa and Lamb (1981) is utilized for the non-linear term computation. If dissipation is necessary, an upwinding scheme is utilized instead. Therefore this study may be interpreted as a machine learning framework for devising hybrid schemes for non-linear term computation with a view to reconstructing turbulence statistics in a superior fashion. Therefore, this study is similar to that employed by Ling and Kurzawski Ling and Kurzawski (2017) for adaptively determining RANS corrections. We note that the classification framework devised in this study is also deployed in an aligned work to switch between functional and structural explicit LES hypotheses spatio-temporally Maulik et al. (2018b) thus proving that high-fidelity DNS statistics may be qualitatively utilized to inform modeling strategies through conditional probability predictions. The article shall describe how the proposed framework is effective in moderating the larger dissipation of an upwinded-scheme through assessments on the Kraichnan turbulence test-case.

II Turbulence modeling equations

The governing equations for two-dimensional turbulence are given by the Navier-Stokes equations in the vorticity-stream function formulation. In this formulation, our non-dimensional governing equation for incompressible flow may be represented as

[TABLE]

where $Re$ is the Reynolds number, $\omega$ and $\psi$ are the vorticity and stream function respectively connected to each other through the Poisson equation given by

[TABLE]

It may be noted that the Poisson equation implicitly ensures a divergence-free flow evolution. The non-linear term (denoted the Jacobian) is given by

[TABLE]

The stream function and the two-dimensional velocity components are related as

[TABLE]

A reduced-order implementation of the aforementioned governing laws (i.e., an LES) is obtained through

[TABLE]

where the overbarred variables are now evolved on a grid with far fewer degrees of freedom. Due to the reduction in supported frequencies, the non-linear Jacobian fails to capture inter-eddy interactions at different wavenumbers. If it is assumed that the finer scales of vorticity are generally dissipative in nature for two-dimensional turbulence (based on Kraichnan’s cascade of enstrophy Kraichnan (1967)), dissipative models may be embedded into the coarse-grained evolution of the vorticity evolution equation to recover some portion of the effect of the finer scales. Explicit LES closures embed dissipation into the vorticity evolution in the form of eddy-viscosity phenomenology or through structural arguments of scale-similarity. However ILES manipulates the computation of the non-linear Jacobian term to add numerical dissipation to mimic that of the unresolved frequencies. The latter framework, while numerically robust, suffers from difficulties associated with directed dissipation where it is often very easy to be over-dissipative in regions where sub-grid dissipation may not be as pronounced. In this article, we introduce a hybrid ILES framework that focuses upwinding at areas where high probability of sub-grid dissipation necessity is detected.

III Non-linear Jacobian computation

The study utilizes two types of non-linear term computation schemes. Our first choice is symmetric, second-order accurate and conserves energy and enstrophy to minimize numerical dissipation. This is given by the well-known second-order Arakawa scheme Arakawa and Lamb (1981) as detailed below. The non-linear term in Equation 5 may be numerically calculated on a coarse grid using

[TABLE]

where $J^{A}(\bar{\omega},\bar{\psi})$ will henceforth refer to the Arakawa discretization. The individual terms on the right hand side of the above equation are given as

[TABLE]

[TABLE]

[TABLE]

The aforementioned scheme is utilized when our proposed classifier recognizes that no dissipation is necessary.

A numerically dissipative computation of the non-linear term allows for that stabilization of noise accumulation at the grid cut-off wavenumbers. Although there are many different methodologies for upwind based dissipation with varying degrees of complexity, in this article, we utilize a conventional upwind-biased scheme as detailed in the following Hoffmann and Chiang (2000). Our ILES Jacobian is computed as

[TABLE]

where

[TABLE]

In addition,

[TABLE]

Note that velocity components are recovered using

[TABLE]

where the second-order accurate reconstruction of the velocity leads to overall second-order accuracy for non-linear Jacobian reconstruction using the upwinded procedure outlined above. We also note that our Poisson equation given by Equation 2 is solved using a spectrally-accurate scheme.

With the choice of one of the two aforementioned schemes, a point in space-time may or may not have an artificial dissipation imparted to it numerically. However, we mention the caveat that switching between these two schemes would mean that the kinetic energy and enstrophy preserving property of the Arakawa scheme is lost.

IV Machine learning for scheme selection

We now discuss the procedure of utilizing DNS data for learning to classify one of the two dissipation scenarios. Of these two options, one is given by the choice of the Arakawa scheme and the other by our upwinded computation of the Jacobian (i.e., when the classification framework has determined that the point does not require sub-grid dissipation or vice-versa respectively). This switching between scenarios is spatio-temporally dynamic. We proceed by outlining our training strategy through the utilization of DNS data. Five equidistant snapshots of DNS data at $Re=32000$ (i.e., at $t=0,1,2,3,4$) and at $N^{2}=2048^{2}$ degrees of freedom (from 40000 available snapshots) are utilized to compute the grid-filtered variables (denoted FDNS) (at $N^{2}=256^{2}$ degrees of freedom) through the application of a spectral cut-off filter. Perfect closure values

[TABLE]

are then obtained (the reader is directed to (Maulik et al., 2019) for details related to the calculation of these quantities). Note here, that the Kraichnan turbulence problem is transient with the evolution of vorticity represented in Figure 1 representing different closure needs over time evolution.

We proceed by introducing the a priori eddy-viscosity given by

[TABLE]

where the right-hand side of the above equation may be calculated from DNS snapshots. The a priori eddy-viscosity is centered at zero (corresponding to where closure modeling is unnecessary) and spreads out in the negative and positive directions (a hallmark of isotropic turbulence). We segregate this a priori estimate of sub-grid effects into three categories as follows. The a priori eddy-viscosities calculated from the DNS data are compared with a Gaussian distribution where values lying less than a distance of 1% of the standard-deviation from the mean (which is zero) are labeled as those requiring no dissipation (due to the low strength of the a priori eddy-viscosity). For posterity, we label these points as $k=1$. Positive values lying beyond this range are labeled as those requiring sub-grid dissipation and are labeled $k=2$. Negative values less than 1% of the standard-deviation are also considered to require no dissipation and are labeled $k=3$. This three-category segregation stems from a learning hypothesis that seeks to identify regions in a flow evolution that require structural, functional or no-closure modeling hypothesis. We link labels of negative or nearly-zero eddy-viscosities to the Arakawa classification and positive eddy-viscosities to the upwinded classification. The positive eddy-viscosity prediction would indicate that the sub-grid term at a point is predominantly dissipative in nature at which point the numerical dissipation of the upwinded scheme would be utilized. We note here that the concept of an a priori eddy-viscosity lies firmly within the explicit LES hypothesis. The classifier is therefore instrumental in moderating ILES deployments through a decision making process that recognizes the dissipative (or forcing) nature of the sub-grid quantities.

We note that the choice of 1% as the decision parameter for switching between hypothesis is motivated by a sensitivity study that showed the highest classification accuracy for the ANN framework. Larger choices of this hyper-parameter would result in a classifier that would be prone to classify most points in the ‘no-model’ zone. However, we clarify that the choice of this value is also correlated with the architecture of the ANN. A potential extension of the proposed hypothesis is to combine architecture search algorithms with varying value of the decision hyper-parameters for larger classification accuracies. In addition, the three-category framework is derived from an aligned study Maulik et al. (2018b) where sub-grid models are determined according to negative, positive and nearly-zero a priori eddy-viscosities and utilizes the same learning. This enables use to determine a unified framework for switching between turbulence model hypotheses as well as numerical dissipation scenarios. However, we would like to emphasize that, for the purpose of switching between the Arakawa and upwinded Jacobian computation, a simple two-class framework would also suffice.

A one-hot labeling of our eddy-viscosity classes is utilized for a classification deployment and a schematic for this hypothesis segregation and labeling is shown in Figure 2. The labels indicate the conditional probability of a point belonging to each possible class. As such, the training labels are given by a value of 1 for the particular class that a point belongs to and zeros for other choices. This is because there is no ambiguity in the class a training sample belongs to. Each label for the a priori eddy-viscosity is also associated with a corresponding input kernel of grid-resolved quantities. This kernel is given by a local stencil of vorticity and stream function. There are 9 inputs each for vorticity and stream function given by a query of the field quantity at a point on the coarse grid, 4 adjacent points in each dimension ($x,y$) and the 4 diagonally adjacent points. Each sample of our training data thus consists of 18 inputs of vorticity and stream function and outputs given by one-hot labels for the choice of closure modeling strategy. We then utilize an ANN to establish a relationship between these inputs and outputs. Mathematically, if our input vector $\mathcal{P}$ resides in a $P$-dimensional space and our desired output $\mathcal{Q}$ resides in a $Q$-dimensional space, this framework establishes a map $\mathbb{M}$ as follows:

[TABLE]

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

[TABLE]

where

[TABLE]

is our input vector for each query of the machine learning framework and where

[TABLE]

is the conditional probability of a Jacobian computation (given by a connection to the explicit closure hypothesis). Note that $i,j$ refer to the spatial indices on the coarse-grid (i.e., the point of deployment). The indices $k=1$ and $k=3$ refer to the Arakawa non-linear Jacobian computation and $k=2$ refers to the upwinded computation instead (see Figure 2). Our optimal map $\mathbb{M}$ is then trained by minimizing the categorical cross-entropy loss-function

[TABLE]

where w are the variable weight and bias parameters of the network, $N$ refers to the total number of samples and $K=3$ is the total number of classification scenarios (i.e., negative, positive or nearly-zero a priori eddy-viscosities). Here, $t_{nk}$ refers to the true label of class $k$ and sample $n$ and $y_{nk}$ refers to a corresponding prediction of the learning framework. One-hot encoding ensures that $t_{nk}$ values are always binary Bishop (2006) and the outputs of the ANN may be interpreted as conditional-probabilities. Our optimal architecture is given by five 40-neuron hidden layers (obtained via grid-search hyper-parameter tuning). All hidden layers utilize ReLU units to impart non-linearity to the layer-wise transformations. For reference, our architecture is trained using the open-source deep learning software Tensorflow and is optimized with the use of ADAM, a popular gradient-descent based optimizer Kingma and Ba (2014). Figure 3 shows the progress to convergence for our framework with our optimally trained network displaying approximately 79% accuracy in classifying points to their correct labels. To summarize this section, we train a deep ANN to estimate probabilities of negative, positive or nearly-zero eddy-viscosities which are utilized to decide the choice of the Jacobian computation. We clarify that the decision to deploy a particular hypothesis is obtained by utilization of the classification scenario which has the highest conditional probability.

V Results

V.1 A posteriori deployment

In this section, we detail the results from an a posteriori deployment of the classification framework (denoted ML henceforth) for the Kraichnan test-case. In the LES evolution of the problem, a considerably coarser grid is used (at $N^{2}=256^{2}$). We remark that the forward deployment of our framework needs to overcome the challenge of numerical errors and is a robust test of the generalizability and robustness of our learning. Our LES results are assessed using angle-averaged kinetic energy spectra and through structure functions of vorticity. In addition, qualitative comparisons are also provided through visual examinations of the vorticity contours. We remark that the LES deployment is performed from $t=0$ to $t=4$ which spans the training regime data obtained from DNS. In what follows we note that DNS refers to a high-fidelity evolution of the governing equations (i.e., at $N^{2}=2048^{2}$ degrees of freedom), UNS refers to results obtained using the Arakawa scheme alone and ILES refers to a simulation where the non-linear Jacobian at all points in space and time are upwinded. Figure 4 shows the a posteriori performance of the proposed framework at $Re=32000$ in terms of energy spectra predictions. The reader may find an exact definition of the kinetic-energy spectra in Maulik and San Maulik and San (2017b). We note that the training data was obtained for the same Reynolds number as well. The prediction of the proposed framework is seen to agree remarkably well with DNS. It is apparent that the switching of schemes using the classifier has obtained an optimal balance between both techniques.

Vorticity contours for LES resolution assessments are shown in Figure 5, where it is apparent that the proposed framework optimally balances the energy-conserving and dissipative natures of the Arakawa and upwinded schemes respectively. This is verified by qualitative examination with FDNS contours obtained by spectrally filtering the DNS snapshot for $Re=32000$ at $t=4$.

A second statistically significant quantity of interest studied in this investigation is the vorticity structure function Grossmann and Mertens (1992) given by

[TABLE]

where the angle-brackets indicate ensemble averaging and $x,y$ indicate a position on the grid with $r$ being a certain distance from this location. Figures 6 and 7 show the structure functions obtained from a posteriori deployments of the UNS, ILES and ML frameworks compared against those obtained from the final time FDNS snapshot. It is clear that the proposed framework balances between UNS and ILES deployments well to recover appropriate trends. We can thus claim that our learning is appropriate for hybrid deployments of dissipative and conservative frameworks for two-dimensional turbulence. Before moving on, we would like to point out to the reader here that the proposed methodology for closure does not require any post-processing prior to deployment in the forward simulation as utilized in several data-driven turbulence modeling studiesBeck, Flad, and Munz (2018); Maulik et al. (2019).

V.2 Validation of learning

In this section, we proceed with a rigorous validation of our learning for deployment in regimes that are not a part of the training data. This is to ensure that the framework has truly learnt a classification based on the underlying physical hypothesis used for data segregation and is not memorizing data. This ensures that our classifier can be used in a more generalizable fashion. Figure 8 shows kinetic energy spectra obtained from the forward deployment of the ML framework for a $Re=64000$ which represents a classification task that the framework has not previously seen (although the physics of the test-case remains similar). As observed, the proposed method performs quite well in this out-of-training data range as well. We note that a similar resolution ($N^{2}=256^{2}$) is utilized for this deployment. In contrast, Figure 9 shows the performance of the ML technique for a reduced resolution of $N^{2}=128^{2}$ but utilizing the same Reynolds number of 32000. The kinetic energy spectra show a successful stabilization of the flow evolution at this reduced resolution although some forcing to the large scales is observed. This suggests that the classification framework may be improved by sampling from different resolutions.

VI Concluding remarks

In this article, we have proposed a neural network based classifier that enables us to take decisions on the choice of non-linear term computation in the LES evolution of the Kraichnan turbulence test-case. The classifier outputs conditional probabilities for the presence (or absence) of eddy-viscosities within three different ranges during deployment and is used to switch between the Arakawa and upwind computation of the non-linear Jacobian for a hybrid upwinded deployment that optimally directs dissipation on the coarse-grained flow field. Our machine learning framework is trained by calculating a priori eddy-viscosities which are projected onto a Gaussian distribution and segregated into three categories. Each category is devised to capture a unique behavior of the underlying sub-grid terms with negative and nearly-zero eddy-viscosity classes signifying absence of sub-grid dissipation. An optimally trained classifier is then utilized to identify if a point requires sub-grid dissipation based on if it is placed in the positive eddy-viscosity category. If so, the upwind Jacobian is calculated for imparting numerical dissipation.

We perform a posteriori assessments on the Kraichnan turbulence test-case through statistical quantities such as the angle-averaged kinetic energy spectra and the vorticity structure functions. It is observed that the proposed framework is successful in balancing the dissipative nature of the upwind scheme and the energy-conserving Arakawa scheme to give excellent agreement with DNS statistics. Validation for out-of-training regimes also indicate that the framework is able to learn the link between grid-resolved quantities at a coarse resolution and the nature of the sub-grid forcing.

Our conclusions therefore point toward the possibility of using classifiers for the unified deployment of numerical schemes with varying dissipation through the decision making process described above. A key strength of our hypothesis stems from the fact that an ILES deployment is moderated by concepts drawn from the explicit LES ideology (i.e., that of an a priori eddy-viscosity). The successful deployment of our method thus points towards the possibility of deploying directed numerical dissipation that preserves the statistics of turbulence without sacrificing the shock-capturing ability of many non-oscillatory schemes. Our future work lies in that particular direction.

Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research under Award Number DE-SC0019290. OS gratefully acknowledges their support. Disclaimer: This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.