A deep learning enabler for non-intrusive reduced order modeling of fluid flows

TL;DR

This paper presents a modular deep neural network framework for non-intrusive reduced order modeling of fluid flows, capable of learning system evolution from data without requiring prior knowledge of governing equations.

Contribution

It introduces a flexible DNN-based approach for data-driven reduced order modeling applicable to various numerical schemes and fluid flow examples.

Findings

High accuracy in predicting system evolution

Robustness across different fluid flow scenarios

Effective non-intrusive model reduction

Abstract

In this paper, we introduce a modular deep neural network (DNN) framework for data-driven reduced order modeling of dynamical systems relevant to fluid flows. We propose various deep neural network architectures which numerically predict evolution of dynamical systems by learning from either using discrete state or slope information of the system. Our approach has been demonstrated using both residual formula and backward difference scheme formulas. However, it can be easily generalized into many different numerical schemes as well. We give a demonstration of our framework for three examples: (i) Kraichnan-Orszag system, an illustrative coupled nonlinear ordinary differential equations, (ii) Lorenz system exhibiting chaotic behavior, and (iii) a non-intrusive model order reduction framework for the two-dimensional Boussinesq equations with a differentially heated cavity flow setup at various Rayleigh numbers. Using only snapshots of state variables at discrete time instances, our data-driven approach can be considered truly non-intrusive, since any prior information about the underlying governing equations is not required for generating the reduced order model. Our \textit{a posteriori} analysis shows that the proposed data-driven approach is remarkably accurate, and can be used as a robust predictive tool for non-intrusive model order reduction of complex fluid flows.

Figures

Tables

Table 1. Table 1: The output of different DNN frameworks learned through training. The trained parameter is then used to update the solution in time, starting from the initial condition.

Neural network framework	Predicted variable	Solution update
DNN-S	$𝒓={𝒚}^{\left(n+1\right)}$	${𝒚}^{\left(n+1\right)}=𝒓$
DNN-R	$𝒓={𝒚}^{\left(n+1\right)}-{𝒚}^{\left(n\right)}$	${𝒚}^{\left(n+1\right)}={𝒚}^{\left(n\right)}+𝒓$
DNN-B	$𝒓=\frac{3{𝒚}^{\left(n+1\right)}-4{𝒚}^{\left(n\right)}+{𝒚}^{\left(n-1\right)}}{2\Delta t}$	${𝒚}^{\left(n+1\right)}=\frac{4}{3}{𝒚}^{\left(n\right)}-\frac{1}{3}{𝒚}^{\left(n-1\right)}+\frac{2}{3}𝒓\Delta t$

Table 2. Table 2: Quantitative assessment of different DNN frameworks with different number of inputs to the neural network for the Kraichnan-Orszag system using the total root mean square error given by Equation  6 .

Framework	$\text{RMSE}\left(p=1\right)$	$\text{RMSE}\left(p=4\right)$
DNN-S	$1.71\times {10}^{-2}$	$1.90\times {10}^{-2}$
DNN-R	$1.02\times {10}^{-3}$	$2.48\times {10}^{-4}$
DNN-B	$5.67\times {10}^{-4}$	$5.97\times {10}^{-4}$

Table 3. Table 3: Quantitative assessment of Galerkin projection and different DNN frameworks for vorticity and temperature modal coefficients for Ra = 3.4 × 10 5 3.4 superscript 10 5 3.4\times 10^{5} and Ra = 9.4 × 10 5 9.4 superscript 10 5 9.4\times 10^{5} using the total root mean square error given by Equation ( 6 ).

Framework	$\text{RMSE}\left(a\right)$	$\text{RMSE}\left(b\right)$
Ra = $3.4\times {10}^5$
ROM-G	$9.8\times {10}^{-1}$	$2.1\times {10}^{-2}$
DNN-S ($p=1$)	$9.9\times {10}^{-2}$	$6.3\times {10}^{-3}$
DNN-R ($p=1$)	$4.8\times {10}^{-1}$	$1.3\times {10}^{-2}$
DNN-B ($p=1$)	$2.9\times {10}^{-2}$	$5.3\times {10}^{-3}$
DNN-S ($p=4$)	$9.1\times {10}^{-2}$	$8.8\times {10}^{-4}$
DNN-R ($p=4$)	$1.2\times {10}^{-1}$	$2.2\times {10}^{-4}$
DNN-B ($p=4$)	$2.0\times {10}^{-1}$	$5.9\times {10}^{-4}$
Ra = $9.4\times {10}^5$
ROM-G	$3.1\times {10}^2$	$1.3\times {10}^0$
DNN-S ($p=1$)	$8.9\times {10}^0$	$1.6\times {10}^{-1}$
DNN-R ($p=1$)	$8.3\times {10}^0$	$1.3\times {10}^{-1}$
DNN-B ($p=1$)	$7.7\times {10}^0$	$1.7\times {10}^{-1}$
DNN-S ($p=4$)	$4.9\times {10}^0$	$1.8\times {10}^{-1}$
DNN-R ($p=4$)	$5.6\times {10}^0$	$1.0\times {10}^{-1}$
DNN-B ($p=4$)	$4.6\times {10}^0$	$1.0\times {10}^{-1}$

Table 4. Table 4: Statistics of Nusselt number computed on the left wall ( x = 0 𝑥 0 x=0 ) for Ra = 3.4 × 10 5 3.4 superscript 10 5 3.4\times 10^{5} and Ra = 9.4 × 10 5 9.4 superscript 10 5 9.4\times 10^{5} . The mean Nusselt number and its standard deviation is computed from instantaneous Nusselt number from time period t = 0 𝑡 0 t=0 to t = 100 𝑡 100 t=100 .

Framework	$\mu$	$\sigma$
Ra = $3.4\times {10}^5$
FOM	-4.5425	$2.46\times {10}^{-3}$
True	-4.5494	$2.48\times {10}^{-3}$
ROM-G	-4.5492	$2.55\times {10}^{-3}$
DNN-S ($p=1$)	-4.5494	$2.48\times {10}^{-3}$
DNN-R ($p=1$)	-4.5494	$2.50\times {10}^{-3}$
DNN-B ($p=1$)	-4.5494	$2.49\times {10}^{-3}$
DNN-S ($p=4$)	-4.5494	$2.48\times {10}^{-3}$
DNN-R ($p=4$)	-4.5494	$2.48\times {10}^{-3}$
DNN-B ($p=4$)	-4.5494	$2.48\times {10}^{-3}$
Ra = $9.4\times {10}^5$
FOM	-5.8411	$4.81\times {10}^{-2}$
True	-5.8603	$2.29\times {10}^{-2}$
ROM-G	-6.2530	$2.43\times {10}^{-1}$
DNN-S ($p=1$)	-5.8572	$3.45\times {10}^{-2}$
DNN-R ($p=1$)	-5.8643	$2.14\times {10}^{-2}$
DNN-B ($p=1$)	-5.8639	$2.19\times {10}^{-2}$
DNN-S ($p=4$)	-5.8634	$2.87\times {10}^{-2}$
DNN-R ($p=4$)	-5.8637	$2.50\times {10}^{-2}$
DNN-B ($p=4$)	-5.8642	$2.29\times {10}^{-2}$

Equations

$$S^l={\mathbf{W}}^l{X}^{l-1},$$

$$X^l=\zeta (S^l+B^l),$$

$$Y={\zeta }_L({\mathbf{W}}^L,B^L,\dots ,{\zeta }_2({\mathbf{W}}^2,B^2,{\zeta }_1({\mathbf{W}}^1,B^1,X))),$$

$$\zeta (\chi )=\text{max}(0,\chi ),$$

$$\frac{d\mathbf{y}}{dt}=\mathbf{F}(\mathbf{y},t),$$

$$\text{RMSE}=\sum _{k=1}^R\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(y_k^{(i)}-{\tilde{y}}_k^{(i)}\right)}^2},$$

$$\frac{d{y}_1}{dt}=y_1{y}_3;\frac{d{y}_2}{dt}=-y_2{y}_3;\frac{d{y}_3}{dt}=-y_1^2+y_2^2,$$

$$\begin{array}{rl}\frac{d{y}_1}{dt}& =\alpha (y_2-y_1),\\ \frac{d{y}_2}{dt}& =y_1(\rho -y_3)-y_2,\\ \frac{d{y}_3}{dt}& =y_1{y}_2-\beta y_3.\end{array}$$

$$\frac{\partial \mathbf{u}(\mathbf{x},t)}{\partial t}=\mathfrak{R}(\mathbf{u}(\mathbf{x},t);\mathbf{x},t),$$

$$\mathbf{u}(\mathbf{x},t)=\sum _{k=1}^R{a}_k^{\mathbf{u}}(t){\varphi }_k^{\mathbf{u}}(\text{x}),$$

$${\int }_{\Omega }{\varphi }_i^{\mathbf{u}}(\text{x}){\varphi }_j^{\mathbf{u}}(\text{x})d\mathbf{x}={\delta }_{ij},$$

$$a_k^{\mathbf{u}}(t_n)={\int }_{\Omega }\mathbf{u}(\mathbf{x},t_n){\varphi }_k^{\mathbf{u}}(\text{x})d\mathbf{x}.$$

$$\mathfrak{M}:\{{\mathbf{a}}^{\mathbf{u}}(t_n),{\mathbf{a}}^{\mathbf{u}}(t_{n-1}),\dots {\mathbf{a}}^{\mathbf{u}}(t_{n-p+1})\}\Rightarrow {\mathbf{a}}^{\mathbf{u}}(t_{n+1}).$$

$$a_k^{\mathbf{u}}(t_0)={\int }_{\Omega }\mathbf{u}(\mathbf{x},t_0){\varphi }_k^{\mathbf{u}}(\text{x})d\mathbf{x}.$$

$$\mathbf{u}(\mathbf{x},t_n)=\sum _{k=1}^R{a}_k^{\mathbf{u}}(t_n){\varphi }_k^{\mathbf{u}}(\text{x}).$$

$$\nabla \cdot \mathbf{u}=0,$$

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\nabla p+\frac{1}{\text{Re}}{\nabla }^2\mathbf{u}+\text{Ri}\theta {\hat{e}}_j,$$

$$\frac{\partial \theta }{\partial t}+(\mathbf{u}\cdot \nabla )\theta =\frac{1}{\text{Re}\text{Pr}}{\nabla }^2\theta ,$$

$$\text{Ra}=\text{Ri}{\text{Re}}^2\text{Pr}.$$

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

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

$$\text{J}(f,g)=\frac{\partial g}{\partial y}\frac{\partial f}{\partial x}-\frac{\partial g}{\partial x}\frac{\partial f}{\partial y}.$$

$$u=\frac{\partial \psi }{\partial y},v=-\frac{\partial \psi }{\partial x}.$$

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

$$\text{Nu($t$)}=\frac{1}{H}\int_{0}^{H}\frac{\partial\theta}{\partial x}\Big{|}_{x=0}dy,$$

$$\bar{\omega }(\text{x})=\frac{1}{M}\sum _{n=1}^M\omega (\text{x},t_n).$$

$${\omega }^{\prime }(\text{x},t_n)=\omega (\text{x},t_n)-\bar{\omega }(\text{x}).$$

$$c_{ij}=⟨{\omega }^{\prime }(\text{x},t_i),{\omega }^{\prime }(\text{x},t_j)⟩,$$

$$⟨f,g⟩={\int }_{\Omega }f(\text{x})g(\text{x})d\text{x}.$$

$$CW=W\mathbf{\Lambda },$$

$${\varphi }_k^{\omega }(\text{x})=\frac{1}{\sqrt{{\lambda }_k^{\omega }}}\sum _{n=1}^M{w}_k^n{\omega }^{\prime }(\text{x},t_n),$$

$${\delta }_{ij}=\{\begin{array}{ll}1,& \text{if }i=j,\\ 0,& \text{if }i\ne j.\end{array}$$

Full text

Dynamical systems, Deep learning, Neural networks, Reduced order modeling, Proper orthogonal decomposition

I Introduction

Many realistic transient flows typically involve very wide ranges of spatial and temporal scales, which place an enormous computational burden on direct numerical simulations (DNS) of such flows based on governing equations. Advancement in high-performance computing systems along with the development of consistent, stable, convergent numerical schemes, and efficient parallel algorithms has enabled us to analyze and study complex real world processes. For instance, it is now possible to collect very high-resolution DNS data relevant to selected turbulent flows which cannot be gathered experimentally Ishihara, Gotoh, and Kaneda (2009). However, the computational cost of performing DNS scales roughly as $\text{Re}^{3}$, where Re is the Reynolds number of the flow Karniadakis and Orszag (1993). Hence, with the present state of the art computing architectures Reed and Dongarra (2015), such high-resolution simulations of multiphysics flows might require weeks of computations even for simple geometries. The situation worsens when a series of numerical simulations need to be run for any parametric design optimization study. To alleviate this, coarse-graining approaches, as performed for example in large eddy simulations (LES), are commonly used to reduce this computational burden Mason (1994); Piomelli (1999); Meneveau and Katz (2000); Moin (2002).

However, the computational cost of full-order simulations (i.e., DNS or even LES) can still be considered extremely prohibitive due to a large number of degrees of freedom needed to resolve all of the flow features, especially in settings where the traditional methods require repeated model evaluations over a large range of parameters. Therefore, many successful model order reduction approaches have been introduced Lucia, Beran, and Silva (2004); Benner, Gugercin, and Willcox (2015); Brunton and Noack (2015); Taira et al. (2017); Rowley and Dawson (2017). The main purpose of such approaches is to reduce this computational burden and serve as surrogate models for efficient computational analysis of fluid systems. A common objective in such reduced order modeling (ROM) approaches is to determine how well these approaches can reproduce the flow dynamics.

Intrusive finite dimensional low order models routinely arise when we apply Galerkin type projection techniques to infinite dimensional modelsLakshmivarahan and Wang (2008a, b); Wang and Lakshmivarahan (2009). On the other hand, without prior information on the governing equations, their operator forms or parameterizations to account for complex physical processes, a non-intrusive ROM approach can be reconstructed to infer such underlying physics from the data itself. Having ability to facilitate dynamic data exchange easier between different components, these non-intrusive models can arguably be more promising and impactful in numerous interdisciplinary fields. Moreover, with the advent of digital twin technologiesTao et al. (2019), the collection of data from sensors has become possible at different stages of product’s lifecycle, and model order reduction might be considered as a key enabler for this digital twin vision in many emerging cyber-physical systemsHartmann, Herz, and Wever (2018).

Reduced order models offer promises in many fields such as system identification Juang and Pappa (1985); Harish and Kumar (2016); Fung and Ruthotto (2019), control Ito and Ravindran (1998); Kunisch and Volkwein (1999); Hovland, Gravdahl, and Willcox (2008); Bergmann and Cordier (2008); Noack, Morzynski, and Tadmor (2011); Kärcher et al. (2018), optimization Lassila and Rozza (2010); Benner, Sachs, and Volkwein (2014); Heinkenschloss and Jando (2018), and data assimilation Cao et al. (2007); Daescu and Navon (2008); Ştefănescu, Sandu, and Navon (2015) applications. In these model reduction approaches, we aim at obtaining simplified (but dense) models from high-fidelity numerical simulation data or data collected from the experiment Rowley and Dawson (2017); Jolliffe and Cadima (2016). To fulfill their objectives in multiple forward simulations of the problem with different model parameters Allen, Weickum, and Maute (2004); Frangos et al. (2010); Degroote, Vierendeels, and Willcox (2010); Weickum, Eldred, and Maute (2006); Shenefelt et al. (2002); Lieberman, Willcox, and Ghattas (2010), these models should be sufficiently accurate and computationally much faster than the high-fidelity numerical simulation. Therefore, there has been progress made in recent years to develop such ROM approaches specifically for nonlinear systems Everson and Sirovich (1995); Chaturantabut and Sorensen (2010); Noack et al. (2003); Carlberg, Bou-Mosleh, and Farhat (2011); Alla and Kutz (2017); Loiseau, Noack, and Brunton (2018).

The basic philosophy of projection-based ROM approaches is to reduce the high degrees of freedom of a governing equation through an expansion in a transformed space, traditionally with orthogonal basis. Among the large variety of projection-based ROM strategies, the proper orthogonal decomposition (POD) has emerged as a popular technique for the study of dynamical systems Berkooz, Holmes, and Lumley (1993); Antoulas, Sorensen, and Gugercin (2000); Benner, Gugercin, and Willcox (2015); Rowley and Dawson (2017), which targets the most dominant characteristics of the flow considering the largest energy containing modes. The POD technique was first introduced in fluid community in the context of extracting coherent structure from turbulent flow field Lumley (1967). Several methods have also been proposed in literature aimed at improving the POD modes Bui-Thanh et al. (2007); Hay, Borggaard, and Pelletier (2009); Kunisch and Volkwein (2010); Carlberg and Farhat (2011). There are also different variants of POD that have been introduced, like spatio-temporal biorthogonal decomposition Aubry, Guyonnet, and Lima (1991), spectral POD (SPOD)Sieber, Paschereit, and Oberleithner (2016), frequency based POD that is also called as SPODTowne, Schmidt, and Colonius (2018), multiscale POD (MPOD)Mendez, Balabane, and Buchlin (2019) which splits the correlation matrix into the contribution of different scales.

The evolution equations for the lower order system are then obtained using the Galerkin projection method. For many flows, the POD-Galerkin method provides an efficient and accurate way to generate ROM methodologies Aubry et al. (1988); Rowley, Colonius, and Murray (2004); Barone et al. (2009); Akhtar, Nayfeh, and Ribbens (2009); Kalb and Deane (2007); Bergmann, Bruneau, and Iollo (2009); Ravindran (2000); Amsallem et al. (2009); San and Borggaard (2015). Furthermore, several successful closure models have been suggested in order to model the effects of discarded modes Borggaard, Iliescu, and Wang (2011); Wang et al. (2012); San and Iliescu (2014); Östh et al. (2014). The POD-Galerkin intrusive approach can also be stabilized with a nonlinear eddy viscosity modelCordier et al. (2013) or with proper selection of linear quadratic coefficientsSchlegel and Noack (2015). However, the projection-based model reduction approaches have limitations especially for complex systems such as general circulation models, since there is a lack of access to the full-order model operators or the complexity of the forward simulation codes that render the need for obtaining the full-order operators Lassila et al. (2014); Peherstorfer and Willcox (2016); Siddiqui et al. (2019).

One of the challenges in Galerkin projection is the deformation of POD modes. As recently discussed by Reiss et al.Reiss et al. (2018), transport-dominated phenomena are usually a challenge for modal methods, since their dynamics cannot be captured accurately by a few dominant spatial modes. If we include more number of modes to better recover the embedded structures in the underlying system, the computational expense increases and the ROM might not be efficient from practical point of view. Furthermore, the construction of a least-order state space is crucial especially in control since every degree of freedom can amplify noiseEhlert et al. (2019). A data-driven manifold learning model has been proposed as a general dynamic ROM modeling frameworkLoiseau, Noack, and Brunton (2018). Ehlert et al.Ehlert et al. (2019) have also presented a manifold representation of the transient oscillatory cylinder wake using a locally linear embedding approach as encoder. They found that this representation outperforms a 50-dimensional POD expansion from the same data. Another key dynamic problem is that hyperbolic convection problems are treated with an elliptic Galerkin methodNoack (2016). Even though the flow has certain specific direction, the Galerkin projection assumes that the modes are globally coupled. This mismatch between Navier-Stokes equations and Galerkin dynamics might not be curable. Also, the frequency range of high-dimensional Navier-Stokes solutions is not resolved in the low-dimensional deterministic systemAubry et al. (1988); Bourgeois, Noack, and Martinuzzi (2013). RempferRempfer (2000) has demonstrated some implications of projecting the Navier-Stokes equations onto low-dimensional bases and showed how the restriction to a low-dimensional basis as well as improper treatment of boundary conditions might affect the validity of ROM. Due to all these limitations, there is a recent interest in generating fully non-intrusive approaches without the need for access to full-order model operators to establish surrogate models Audouze, De Vuyst, and Nair (2013); Mignolet et al. (2013); Xiao et al. (2015, 2015); Peherstorfer and Willcox (2016); Hesthaven and Ubbiali (2018); Hampton et al. (2018); Chen et al. (2018a); Xiao et al. (2019); Wang, Hesthaven, and Ray (2019). Dynamic mode decomposition (DMD) modelsRowley et al. (2009); Schmid (2010); Jovanović, Schmid, and Nichols (2014); Hemati, Williams, and Rowley (2014); Proctor and Eckhoff (2015); Noack et al. (2016) provide this non-intrusive representation directly from data by their nature, and several approaches have been readily available for optimal mode selectionBistrian and Navon (2015); Alekseev et al. (2016); Bistrian and Navon (2017); Pascarella, Fossati, and Barrenechea (2019).

There is a broad range of opportunities for the application of machine learning algorithms to develop non-intrusive reduced order models. A good discussion on the application of data-driven methods for dynamical systems can be found in a book by Brunton and KutzBrunton and Kutz (2019). We also refer to a recent review article Brunton, Noack, and Koumoutsakos (2019) for a comprehensive overview of the machine learning literature in fluid mechanics. A number of studies have been done to apply data-driven techniques to predict the high-dimensional complex dynamical systems San and Maulik (2018a); Xie, Zhang, and Webster (2018); Raissi, Perdikaris, and Karniadakis (2018); Chen et al. (2018b); Rudy et al. (2019); Kani and Elsheikh (2017); Yeo and Melnyk (2019); Qin, Wu, and Xiu (2019); Rahman, Rasheed, and San (2018). San and MaulikSan and Maulik (2018b) proposed a methodology to account for the effects of truncated POD modes using a single layer feed-forward neural network. A multistep neural network was proposed to identify nonlinear dynamical system from the data by combining classical numerical analysis techniques with the powerful nonlinear approximation capability of neural networks Raissi, Perdikaris, and Karniadakis (2018). Xie, Zhang, and WebsterXie, Zhang, and Webster (2018) used the multistep neural network to approximate the full order model projected on low-dimensional space with a supervised learning task. A deep residual recurrent neural network was introduced as an efficient model reduction technique for nonlinear dynamical systems Kani and Elsheikh (2017). Vlachas et al.Vlachas et al. (2018) developed the data-driven forecasting method for high-dimensional, chaotic systems using the hybrid approach which combines mean stochastic model and recurrent long short-term memory (LSTM) neural network. The LSTM recurrent neural network was used to model the temporal dynamics of turbulence in a ROM framework Mohan and Gaitonde (2018). Pathak et al.Pathak et al. (2018) proposed a hybrid forecasting model combining the knowledge of the governing equation of the dynamical system and machine learning technique to predict the long term behavior of chaotic systems. This hybrid approach was found to be better than either its pure data-driven component or its model-based component.

In our proposed ROM framework, we will bypass the Galerkin projection step of the projection based ROM with our proposed neural network architectures to build a fully non-intrusive approach. This non-intrusive ROM (NIROM) framework can be viewed as a decomposition of the problem into basis representation and forecasting subproblems. We illustrate our NIROM approach using the deep feed-forward neural network architectures. However, it can be easily applied to other types of neural networks (as demonstrated for data-driven forecasting of dynamical systems Vlachas et al. (2018); Lu et al. (2018)) or more traditional time series forecasting tools Shumway and Stoffer (2017). The neural networks are capable of approximating the nonlinear functions and have been successfully used in turbulence modelingLing, Kurzawski, and Templeton (2016); Milano and Koumoutsakos (2002); Gamahara and Hattori (2017), solving differential equationRaissi, Perdikaris, and Karniadakis (2019); Michoski et al. (2019). We learn the dynamics of the reduced order model directly from the output of the full order model projected on the low-dimension space using a supervised learning task. The main advantage of this non-intrusive approach is that it does not require the information about the equations governing the full order model. Although the proposed approach helps to generate a NIROM framework solely from the snapshot data reconstructed onto a POD-spanned space, it may still suffer from fundamental challenges of traditional POD-Galerkin models (e.g., we refer to Zerfas et al.Zerfas et al. (2019) for a recent discussion about ways to mitigate their lack of accuracy).

The paper is organized as follows: Section II introduces deep neural network architecture and implementation of different DNN frameworks for dynamical systems. Section III gives numerical results using our DNN frameworks for two dynamical systems: Kraichnan-Orszag system and Lorenz system. We present the generalized non-intrusive ROM framework in Section IV.1. In Section IV.2, we present the Boussinesq equation problem and its POD analysis. The non-intrusive ROM framework for the Boussinesq equation is described in Section IV.3. The numerical results in Section V demonstrate the effectiveness of our non-intrusive approach for reduced order modeling of differentially heated cavity problem at two Rayleigh numbers. We give some concluding remarks and suggestions for future work in Section VI.

II Learning Framework

The deep neural network is an artificial neural network composed of several layers made up of the predefined number of nodes. These nodes are also called neurons. A node combines the input from the data with a set of coefficients called weights. These weights either amplify or dampen the input and thereby assign the significance to the input with respect to the output that the DNN is trying to learn. In addition to the weights, these nodes have a bias for each input to the node. The input-weight product and the bias are summed and the sum is passed through a node’s activation function. The above process can be described using the matrix operation as given by Hagan et al. (1996)

[TABLE]

where $X^{l-1}$ is the output of the $(l-1)$th layer, $\mathbf{W}^{l}$ is the matrix of weights for the $l$th layer. The output of the $l$th layer is given by

[TABLE]

where $B^{l}$ is the vector of biasing parameters for the $l$th layer and $\zeta$ is the activation function. If there are $L$ layers between the input and the output, then the mapping of the input to the output can be derived as follow

[TABLE]

where $X$ and $Y$ are the input and output of the deep neural network, respectively.

The input layer usually takes the raw data from the training dataset and transfers it to the second layer. It does not have any biasing or activation through an activation function. The output layer usually has the linear activation function and some bias associated with the inputs. The linear activation function simply takes the summation of inputs received from the previous hidden layer and the associated bias of the output layer. In this study, we use ReLU activation function for all hidden layers and linear activation function for the output layer in all DNN frameworks. The ReLU activation function can be expressed as

[TABLE]

where $\zeta$ is an activation function and $\chi$ is the input to the node.

Each entry of the matrix $\mathbf{W}$ and $B$ is learned through backpropagation and some optimization algorithm. The backpropagation algorithm provides a way to compute the gradient of the objective function efficiently and the optimization algorithm gives a rapid way to learn optimal weights. The objective of the neural networks in this study is to learn the weights associated with each node in such a way that the root mean square error between the true labels $Y_{0}$ and output of the neural network $Y$ is minimized. The backpropagation proceeds as follows: (i) the input and output of the neural network is specified along with initial weights, (ii) the training data is run through the network to produce output $Y$ whose true value is $Y_{0}$, (iii) the derivative of the objective function with each of the training weight is computed using the chain rule, (iv) the weights are updated based on the learning rate and then we go to step (ii). We continue to iterate through this procedure until convergence or the maximum number of iterations is reached. The Adams optimization algorithm Shahriari et al. (2015) is used in this study for learning optimal weights to minimize the objective function.

Deep neural networks are capable of approximating nonlinear dynamical systems as shown in many studiesRaissi, Perdikaris, and Karniadakis (2018); Yeo and Melnyk (2019); Chen et al. (2018b); Haber and Ruthotto (2017). The general nonlinear dynamical system can be presented by an equation of the form

[TABLE]

where $\bm{y}(t)$ is the state variable at time $t$ and $\bm{F}(\bm{y},t)$ is the nonlinear function evaluated for each component of state variable $\bm{y}(t)$.

Figure 1 shows three different DNN frameworks used in this study to predict the dynamical system. Although their characteristics on the stability and well-posedness are beyond the scope of the present work, we refer to the study of Chang et al.Chang et al. (2018) for several deep neural networks and their stability issues. Our motivation for different choices of DNNs comes from recent development of approximation and discovery of dynamical systems using deep learning techniques Qin, Wu, and Xiu (2019); Raissi, Perdikaris, and Karniadakis (2018). At this point, it may be inconclusive to say which DNN framework is superior as compared to others. However, our empirical evidences suggest that learning residual information or numerical slope helps in more accurate prediction of dynamical systems. In case of the DNN-S framework, we train the neural network to learn an update formula which advances the state of the system from $\bm{y}^{(n)}$ to $\bm{y}^{(n+1)}$ directly, where $(n)$ denotes the state of the system at time $t_{n}$. The S in the DNN-S stands for sequential. The past history of the state of the system can be utilized to predict the system’s future state by incorporating it in the input features of the neural network. If we want to include the state history for $p$ time steps, then the input of the neural network consists of $\bm{y}^{(n)}$, $\bm{y}^{(n-1)},\dots,\bm{y}^{(n-p+1)}$. Hence, the neural network will have $R\times p$ input features and $R$ output labels (i.e., $\mathfrak{M}:{\rm I\!R}^{R\times p}\Rightarrow{\rm I\!R}^{R}$), where $\mathfrak{M}$ refers to the DNN model, and $R$ is the number of components of the dynamical system. Once the neural network is trained and the weights are learned, the neural network is used to predict the state of the system starting with an initial condition $\bm{y}^{(0)}$ and proceeding in time iteratively. The prediction of the system in iterative fashion is shown in Figure 1a. If the neural network is trained using the state of the system for $p$ time steps, then the previous $p$ time steps should be stored in the prediction routine. The future state of the system $\bm{y}^{(n+1)}$ is predicted using the true state of the system for first $p$ time steps. After $p$ time steps, only the predicted values by the neural network are used in the input. The predicted variable of the neural network and the solution update formula during an iterative prediction are given in Table 1.

For the DNN-R framework, we learn the difference between the state of the system at time step $t_{n}$ and next time step $t_{n+1}$. The residual between two time steps is then applied to update the current state during prediction. The learning of the residual information instead of sequential update formula helps in stabilizing the neural networkHaber and Ruthotto (2017) and also improves the accuracy of neural network predictionLu et al. (2018). The DNN-R framework can also be implemented using the history of the system’s state, similar to the DNN-S framework. The related framework was employed for the model order reduction of parametric viscous Burgers equation problemSan, Maulik, and Ahmed (2019) with one temporal leg history (i.e., $p=1$). They call it POD-ANN-RN framework and this framework was found to give better results for interpolatory and extrapolatory ROM than sequential learning. The iterative prediction of the dynamical system using DNN-R framework is shown in Figure 1b. The future state of the system is computed using the solution update formula mentioned in Table 1.

We also introduce an additional framework called a DNN-B framework as shown in Figure 1c. The B in the name stands for the backward difference. In this framework, we use the second order backward difference numerical scheme to compute the slope at time step $t_{n}$. The numerical slope information is then used to update the state of the system during prediction. The similar approach is also implemented in other data-driven methods for dynamical systems. One such work is the Multistep neural network Raissi, Perdikaris, and Karniadakis (2018) used for the data-driven discovery of nonlinear dynamical systems from experimental measurements of the state of the system. In their work, the evolution of the system $\bm{F}$ is learned using the neural network by incorporating the multistep Adams-Moulton numerical scheme in computing the loss function of the neural network. In our work, we use the numerical slope as the predicted variable and use mean squared error as the loss function. One of the advantages of directly learning the numerical slope is that standard loss functions available in Keras library can be directly applied without any modification. We apply the second order backward difference scheme in DNN-B framework to compute the discrete numerical slope. However, the framework can be implemented with any family of numerical schemes such as central difference or forward difference family. The equation used to determine the numerical slope and the solution update formula during iterative prediction are provided in Table 1.

The quantitative performance of each DNN framework is measured by a quantity root mean square error (RMSE). The root mean square is determined for each component $k$ between true state and the state predicted by the neural network from the initial time to final time. The root mean square of each component is added to get the total root mean square error. The root mean square error is defined as:

[TABLE]

where $R$ is the total number of components of the dynamical system, $N$ is the total number of time steps in the evolution of the dynamical system, $\bm{y}$ is the true solution and $\tilde{\bm{y}}$ is the solution predicted by the neural network.

III Time Series Prediction for Dynamical Systems

Before implementing our proposed DNN frameworks within the non-intrusive ROM setup, we demonstrate the capability of our DNN frameworks to model nonlinear dynamical systems using two examples. Section III.1 provides numerical results for three mode Kraichnan-Orszag problem and Section III.2 gives results for the chaotic Lorenz system.

III.1 Kraichnan-Orszag System

We start by considering the nonlinear dynamical system Kraichnan-Orszag Wan and Karniadakis (2005) of order $R=3$ as the first test problem. The Kraichnan-Orszag system is defined as:

[TABLE]

with an initial condition $y_{1}(0)=1$, $y_{2}(0)=0.1\xi$, and $y_{3}(0)=0$. The input $\xi$ lies between the interval $[-1,1]$. The modeling of this system is particularly challenging due to the discontinuity at planes $y_{1}(0)=0$ and $y_{2}(0)=0$.

The system is solved numerically using the SciPy function odeint for time integration between the time interval $[0,10]$ and $N$ = 1000 time steps ($\Delta t=0.01$). The initial condition considered for the system is $\xi=0.5$ (i.e., $[y_{1}\;y_{2}\;y_{3}]^{T}$ = $[1.0\;0.05\;0.0]^{T}$). The true solution is used for training the neural network. We apply same hyperparameters for all DNN frameworks. We use three hidden layers with 128 neurons each. The maximum number of iterations is set to 600 and 10% of the data is used for validation to avoid overfitting.

The true state of the system and the state predicted by different DNN frameworks are shown in Figure 2 for $p=1$. All DNN frameworks are correctly able to predict all three states of the dynamical system. We also see that DNN-R and DNN-B frameworks perform better than the DNN-S framework. For DNN-S framework, the state predicted by the neural network is very slightly shifted from the original state. Figure 3 provides the numerical results for DNN frameworks for $p=4$. The state predicted by all DNN frameworks is almost the same as the true state. Table 2 compares the RMSE calculated using Equation (6) for all three DNN frameworks with different number of inputs to the neural network.

III.2 Lorenz System

The Lorenz systemLorenz (1963) can be described by equations given below

[TABLE]

For the Lorenz system, we use $\alpha=10$, $\rho=28$, and $\beta=8/3$. The modeling of the dynamics of the Lorenz system is a challenging problem due to its highly nonlinear and chaotic behavior. This system arises in many simplified models for physical processes Berg, Pomeau, and Vidal (1986); Poland (1993).

The training data for the neural network is generated using the true solution to the Lorenz system. We use the initial condition $[y_{1}\;y_{2}\;y_{3}]^{T}=[-8\;7\;27]^{T}$ and generate the true solution numerically using the SciPy function odeint. The true solution is generated between the time interval $[0,25]$ with the time step $\Delta t=0.01$ ($N=2500$ time steps). We apply the same hyperparameters for all DNN frameworks. We use five hidden layers with 128 neurons each. The maximum number of iterations is set to 1000 and 10% of the training data is used for validation to avoid overfitting.

Figure 4 shows the true trajectory of the Lorenz system and the trajectory predicted by different DNN frameworks using $p=1$ in the input training data. The time period for which the predicted trajectory is the same as the true trajectory varies for different DNN frameworks. Lorenz system has a chaotic behavior and a smaller error in the predicted state of the system can lead to larger error in the forecasted state of the system. The time period for which the predicted trajectory follows true trajectory is longer for the DNN-R and DNN-B frameworks than the DNN-S framework. Figure 5 shows the similar results for $p=4$ in input training data. If we compare Figure 4a and Figure 5a we see that the predicted trajectory follows the true trajectory longer as we include the temporal history of the state of the system in the input to the neural network. We do not see a similar behavior for DNN-R and DNN-B frameworks.

Even though the neural network is not able to predict the true trajectory of the individual state after some time, we should check the ability of the neural network to capture the overall dynamics of the Lorenz attractor. This is due to the fact that the Lorenz system has a positive Lyapunov exponent and a small perturbation in the initial condition can cause the system to diverge exponentially. Hence, a number of studies checks for the dynamics of the Lorenz attractor rather than the individual state of the system Raissi, Perdikaris, and Karniadakis (2018); Brunton, Proctor, and Kutz (2016). Figures 6 and 7 show the exact dynamics of the Lorenz attractor and predicted dynamics for different DNN frameworks with $p=1$ and $p=4$, respectively. All DNN frameworks are capable of capturing the dynamics of the Lorenz system and accurately predicts the correct shape of the Lorenz attractor.

IV Non-intrusive Reduced Order Modeling (NIROM)

We observe a successful predictive performance for a simple nonlinear dynamical system (governed by 3 coupled ordinary differential equations) using different DNN frameworks in the previous Section. We also saw the ability of the neural network to capture the overall dynamics of the chaotic nonlinear Lorenz system. This has motivated us to test the proposed DNN frameworks for the model order reduction of a real-world test problem.

Indeed, we transform a partial differential equation system into a set of ordinary differential equations in order to get the reduced order model. Our main motivation of this study is to convey the message that the neural network architectures can be used in developing a robust and efficient non-intrusive reduced order model. With a goal to recover the reduced order model dynamics of the underlying flow phenomena, in this Section, we implement all DNN frameworks introduced in Section II in model order reduction of a differentially heated cavity problem. This test problem setup is well-established as a model validation test case due to the problem’s simplicity in terms of problem definition as well as the wide variety of applications such as nuclear reactor core isolation, solar energy storage and so on Paolucci (1994); Podvin and Le Quéré (2001); Johnston and Krasny (2002). At first, we will describe our non-intrusive ROM framework. Then, we will define our test problem setup along with the governing equations, then we will briefly demonstrate our non-intrusive ROM framework for the Boussinesq equation problem. Finally, we will demonstrate the comparative performance of the aforementioned DNN frameworks in terms of the temporal evolution of vorticity and temperature field and contours of temperature field for a given initial condition in Section V.

IV.1 NIROM framework

To develop our NIROM framework, we define a generalized PDE system as below

[TABLE]

where $\mathbf{u}$ refers to the problems of interest (e.g., velocity, pressure, and temperature, etc) and $\mathfrak{R}$ converts the physical process possibly with linear, nonlinear, and forcing terms.

In our NIROM framework, we assume that we do not have access to $\mathfrak{R}(\mathbf{u}(\mathbf{x},t);\mathbf{x},t)$ operator. We formulate the proposed NIROM framework assuming that we have discrete snapshots of $\mathbf{u}(\mathbf{x},t)$. The details steps of our NIROM framework are outlined in Algorithm 1. We highlight that this physics-agnostic modeling approach is quite modular and decomposes the problem into the basis representation and forecasting problems. Figure 8 depicts various stages of NIROM framework for any generalized PDE system.

IV.2 Problem definition: Boussinesq equations

For this numerical experiment, we investigate the convective flow behavior in a two-dimensional buoyancy-driven flow in a differentially heated tall cavity Liu, Wang, and Johnston (2003). We consider the following dimensionless form of the two-dimensional incompressible Boussinesq equations San and Borggaard (2015); Deane and Sirovich (1991); Sirovich and Deane (1991); San and Maulik (2018c) on our computational domain, $\Omega$:

[TABLE]

where $\bf u$ refers to the velocity vector, $p$ and $\theta$ denote the pressure and temperature fields, respectively, and $\hat{e}_{j}$ is the unit vector in $y$ direction. Here, $\nabla$ and $\nabla^{2}$ are the standard two-dimensional differential and Laplacian operator, respectively. In general, the Boussinesq equations are characterized by three dimensionless numbers: Reynolds number (Re), Prandtl number (Pr) and Richardson number (Ri). However, other relevant dimensionless numbers can be introduced into the equation as a control parameter based on the physics of the system. For example, we introduce Rayleigh number (Ra) in our study due to the natural convection heat transfer which can be expressed as

[TABLE]

In our problem setup, we fix Pr$=0.71$, Ri $=1$. In our two-dimensional full order model (FOM) simulation, we utilize the vorticity-streamfunction formulation to avoid numerical complexity associated with the primitive variable formulation San and Maulik (2018c). Therefore, we introduce the vorticity ($\omega=\nabla\times\bf u$) and streamfunction ($\psi$) for Equation (16) specified by the following coupled equations Johnston and Krasny (2002); San and Borggaard (2015):

[TABLE]

where Jacobian J accounts for the nonlinear advection term, which is defined as

[TABLE]

The flow velocity components can be found from the stream function, $\psi$, using following definitions:

[TABLE]

The kinematic equation connecting the vorticity and streamfunction can be found by substituting the velocity components in terms of stream function, which form the following divergence-free constraint satisfying the Poisson equation:

[TABLE]

Our Cartesian computational domain is $(x,y)\in[0,\ 1]\times[0,\ 8]$. We utilize wall boundary conditions for all four sides of our computational domain with an adiabatic condition on the top and bottom wall. We imply Dirichlet conditions on left ($\theta=0.5$) and right ($\theta=-0.5$) walls. At the cavity walls, we enforce no-slip boundary conditions by setting zero values for streamfunction and vorticity values calculated by the Briley’s formula Maulik and San (2017). A detailed derivation and discussion on our test problem setup along with the numerical schemes for generating snapshots through FOM simulation can be found in a recent study conducted by San and Maulik San and Maulik (2018c). Note that, even though the simulation is performed for a maximum time of $t=1000$, we perform all our statistical analysis from $t=900$ (after an initial transient period) at $128\times 1024$ grid resolutions, and collect snapshot data only between $t=900$ and $t=950$, and perform our quantitative analyses for the prediction of both in-sample data zone (between $t=900$ and $t=950$) as well as the out-of-sample data zone (between $t=950$ and $t=1000$). For the rest of this paper, we shall refer to this initial state (i.e., time $t=900$ in our physical simulation) as $t=0$ for prescribing initial conditions for ROMs. Therefore, our in-sample data zone spans between $t=0$ and $t=50$, and our out-of-sample data zone spans between $t=50$ and $t=100$. We use the data in the in-sample-zone for training the neural network.

The DNS technique used for performing the FOM simulation is validated by recording the simulation statistics at a designed probe point and comparing it with the studies in the literature Tyliszczak (2016). We simulate the differentially heated cavity problem for four different Rayleigh numbers. The flow is smooth and orderly for lower Rayleigh numbers. As we increase the Rayleigh number, the flow starts getting chaotic and turbulent. We calculate the Nusselt number along the vertical left wall ($x=0$) using the below formula

[TABLE]

where $H=8$. Figure 9 shows the statistics of Nusselt number along the left wall and the designated probe temperature ($x=0.125$ and $y=7.0$). It can be seen that the flow behavior is periodic for both Nusselt number and the temperature history at the probed location at lower Rayleigh numbers. As the Rayleigh number increases, we see that the Nusselt number and probed temperature variation is not periodic due to the turbulent nature of flow. We test our non-intrusive framework for Ra = $3.4\times 10^{5}$ where the flow is periodic and for Ra = $9.4\times 10^{5}$ where the flow is chaotic.

IV.3 NIROM framework for Boussinesq equations

In this Section, we present our non-intrusive ROM setup for the unsteady, incompressible Boussinesq equations given by Equation (20) and Equation (21). In our ROM settings, we first compute the desired set of orthogonal spatial basis functions from a stored high-fidelity data snapshots using proper orthogonal decomposition (POD). We obtain the data snapshots from a high-resolution FOM simulation. Using the pre-computed basis functions, we develop the non-intrusive framework in an encoder-decoder approach using different DNN frameworks proposed in Section II. To compare the predictive performance of the non-intrusive ROMs with respect to the standard intrusive ROM, we develop our intrusive ROM framework (ROM-G) using the Galerkin projection to derive the dynamical model for the POD coefficients Rowley, Colonius, and Murray (2004); Lucia, Beran, and Silva (2004); Iollo, Lanteri, and Désidéri (2000). The implementation of Galerkin projection for the underlying test problem is detailed in the recent work by San and Maulik San and Maulik (2018c) and hence, is not discussed in the present work. Here, we demonstrate the development of the non-intrusive ROM methodologies briefly before proceeding to the numerical results for analyses.

The POD bases for the vorticity field can be constructed from the field variable $\omega(x,y)$ at different time steps which we denote as snapshots, i.e., for M number of snapshots, $\omega(\textbf{x},t_{n})$ are the stored snapshots for $n=1,2,\dots,M$. The time-averaged field can be computed as

[TABLE]

To map the snapshot data to its origin, we then compute the mean-subtracted snapshots or the fluctuating fields by

[TABLE]

This subtraction guarantees that ROM solution would satisfy the same boundary conditions as the full order model Chen, Tu, and Rowley (2012). To simplify the eigenvalue problem necessary for POD bases calculation, we utilize the standard method of snapshots proposed by Sirovich Sirovich (1987) which reduces the larger dimension problem to a much smaller dimension problem. To do so, we construct a $M\times M$ correlation matrix of the fluctuating part ${C}=[c_{ij}]$ which is computed from the inner product of the mean-subtracted snapshots

[TABLE]

where $i$ and $j$ refer to the snapshot indices. The definition of the inner product of any arbitrary two fields $f$ and $g$ can be expressed as

[TABLE]

In the current study, we use the well-known Simpson’s integration rule for a numerical computation of the inner products. Next, an eigendecomposition of $C$ matrix is performed by solving

[TABLE]

where $\mathbf{\Lambda}$ is a diagonal matrix whose entries are the eigenvalues $\lambda_{k}^{\omega}$ of $C$, and $W$ is a matrix whose columns $w_{k}$ are the corresponding eigenvectors. This has been shown in detail in various POD literature (see, e.g., Sirovich (1987); Ravindran (2000)). It should be noted that eigenvalues need to be arranged in a descending order (i.e., $\lambda_{1}^{\omega}\geq\lambda_{2}^{\omega}\geq\dots\geq\lambda_{M}^{\omega}$), for proper selection of the POD modes. The POD modes of vorticity field $\phi_{k}^{\omega}$ are then computed as

[TABLE]

where $w^{n}_{k}$ is the $n$th component of the eigenvector $W$. The scaling factor, $\left(\dfrac{1}{\sqrt{\lambda_{k}^{\omega}}}\right)$ is to guarantee the orthonormality of POD modes i.e., $\langle\phi_{i}^{\omega},\phi_{j}^{\omega}\rangle=\delta_{ij}$. Here, $\delta_{ij}$ is the Kronecker delta defined by

[TABLE]

Similarly, we can compute the POD bases for the temperature field which is $\phi_{k}^{\theta}(\textbf{x})$.

Figure 10 shows the decay of eigenvalues of the correlation matrix for vorticity and temperature field. We can observe that the first 10 modes are able to capture more than 99% of the total energy for both vorticity and temperature at lower Rayleigh numbers (Ra = $3.4\times 10^{5}$ and Ra = $5.4\times 10^{5}$). Also, there is a faster decay of eigenvalues for lower Rayleigh numbers than for higher Rayleigh numbers. For higher Rayleigh numbers, the first 10 modes are not sufficient to capture a major portion of the energy. For example, the first 10 modes of vorticity field capture only 76% and 69% of the total energy for vorticity at Rayleigh numbers Ra = $7.4\times 10^{5}$ and Ra = $9.4\times 10^{5}$, respectively. In order to capture more than 95% of the energy, we will need to consider more than 40 modes. If we include 40 modes, then the Galerkin projection becomes computationally expensive and we lose the benefit of ROM framework. If we include only 10 modes, then the Galerkin projection is unbounded, as we will see in Section V, and it gives a physically wrong solution. We demonstrate that our non-intrusive framework give sufficiently accurate results comparable to the true projection of the FOM solution on reduced order space even with less number of POD modes. It is evident that if we include a higher number of modes to capture the total energy, the true projection of FOM solution on lower dimensional bases will approximate the FOM solution. However, the computational burden will also go up with an increased number of modes. In Figure 11, we provide the contour plots for a few POD basis for the temperature field $\theta$ to indicate the structure of the solution at Ra = $3.4\times 10^{5}$. We can observe that the solution is smooth and periodic for lower Rayleigh number case and only small structures are truncated after 10 POD modes. On the other hand, we observe from Figure 12, that there are still some of the large structures remaining in the 10th mode for higher Rayleigh number case. This means that some of the important flow features are truncated due to consideration of only the first 10 dominant POD modes.

To obtain the non-intrusive ROM, we utilize an encoder-decoder approach to transfer data from full order space to reduced order space and vice versa. During the encoder stage, we transform data from the full order space to the reduced order space by using the following projection for both field parameters:

[TABLE]

So, at initial time $t=0$, we can compute the initial conditions by

[TABLE]

where $\omega(\textbf{x},0)$ and $\theta(\textbf{x},0)$ is the vorticity and temperature field, respectively, specified at initial time. Here, $a_{k}(t)$ and $b_{k}(t)$ are the time-dependent modal coefficients for the vorticity and temperature, respectively. This will form an initial value two equations ODE system similar to the problem in Section III where we require to solve $da_{k}/dt$ and $db_{k}/dt$ until the final time. As discussed earlier, we utilize the different DNN frameworks and time histories as input to predict the temporal evolution of $a_{k}(t)$ and $b_{k}(t)$. For the prediction with sequential time history data, we predict the next time step using the DNN-S framework with one ($p=1$) and four ($p=4$) legs temporal history in the input. We also employ DNN-R and DNN-B frameworks which predicts the residual and the numerical slope, respectively based on one ($p=1$) and four ($p=4$) legs temporal history in the input to the neural network. In the decoder stage, we reconstruct the reduced order solution to the full order solution by using the following definition:

[TABLE]

where $R\ (<<M)$ is the retained most energetic POD modes. Since we do not use any physical equations for the time integration of the solution field, we can say our ROM setup using DNN is fully non-intrusive.

V Numerical results for Boussinesq equations

To evaluate the performance of our DNN frameworks within the underlying non-intrusive ROM setup, we present the time series evolution for the vorticity transport equation. In addition, we compare the temperature field at the final time (i.e., $t=100$) predicted using non-intrusive ROM framework, with FOM simulation and its projection on reduced order space (True projection). We further evaluate the performance of DNN frameworks in predicting engineering quantities of interest such as time-averaged Nusselt number. We also perform the quantitative assessments of different model’s predictive performance in terms of the RMSE calculated using Equation (6). We present our analysis for two Rayleigh numbers: Ra = $3.4\times 10^{5}$ where the flow is smooth and periodic, and Ra = $9.4\times 10^{5}$ where the flow is turbulent and chaotic.

For all DNN frameworks, we use six hidden layers with 120 neurons each. The maximum number of iterations is set to 900 and 10% of the data is used for validation to avoid overfitting. The root mean square error used for evaluating the quantitative performance of DNN frameworks measures the difference between the true data and the predicted data for each mode. Hence, we present the true and predicted trajectories by neural network for only two modes $a_{1}$ and $a_{7}$ for conciseness. Additionally, we compare the modal coefficient trajectory with the POD Galerkin projection (GP) trajectories (i.e., time dependent amplitude coefficients of ROM-G). It should be noted that our criteria for selection of hyperparameters are not necessarily optimal but are based on heuristics (involving different activation functions, number of layers/neurons/iterations etc.) that enable our neural networks to accurately predict time evolution of dynamical systems. We also highlight that there are statistical methods available to select hyperparameters such as Bayesian optimizationSnoek, Larochelle, and Adams (2012).

As illustrated in Figure 13 and for the rest of our analysis in this Section, we display the true solution as a black solid line, and regular POD-Galerkin projection based ROM, i.e., ROM-G solution as a blue dash-dot line. The solution predicted by the neural network is shown by the dashed red line. The training data for the neural network is taken from the time series of modal coefficients obtained by projecting the FOM solution on the reduced order space between $t=0$ to $t=50$. This is consistent with the snapshot data used for generating the POD bases. The training data for the neural network is highlighted using the light gray color in all time series plots. When we test the neural network, we start with an initial condition at $t=0$ and proceed in an iterative fashion as discussed in Section II. Therefore the data between $t=0$ to $t=50$ is in-sample data and $t=50$ to $t=100$ is the out-of-sample data. The neural network has seen the in-sample data during training and hence is expected to give a good prediction for that time period. The question to ask is how does a neural network predict for the out-of-sample data.

Figure 13 shows the time evolution of two modal coefficients of vorticity transport equation for all DNN frameworks using $p=1$ in the input training data for Ra = $3.4\times 10^{5}$. It can be clearly seen that there is a phase difference between the true projection and Galerkin projection for the first modal coefficient $a_{1}$. Also, the Galerkin projection predicts slight amplification in the amplitude of $a_{1}$, especially near the final time. The similar amplification in magnitude is also observed for $a_{7}$ with the Galerkin projection. The prediction by all DNN frameworks is better than the Galerkin projection and the prediction is very close to the true projection of the FOM solution on reduced order space. For the DNN-R framework, there is a slight phase shift in DNN prediction with respect to the true projection as the time proceeds. The prediction for the DNN-R framework can be improved by including the past history of the modal coefficients in the input training data. Figure 14 shows the DNN prediction with $p=4$ in the input training data. We can see that the prediction is improved for the DNN-R framework using $p=4$. The prediction for DNN-S and DNN-B frameworks was already good using $p=1$ and remains the same with $p=4$ also. We see a similar prediction for other modal coefficients also.

The results presented in Figure 13 and Figure 14 were for Ra = $3.4\times 10^{5}$. At lower Rayleigh number, the flow is smooth and orderly, and hence the evolution of modal coefficients was periodic. This is a simple problem with stationary time series and can be solved using simple methods like extreme learning machineHuang, Zhu, and Siew (2006); San and Maulik (2018a). The DNN will be beneficial for higher Rayleigh number case after the onset of turbulence. Figures 15 and 16 show the similar results for Ra = $9.4\times 10^{5}$. The modal coefficients are not periodic due to the chaotic and turbulent nature of flow taking place in the cavity at such higher Rayleigh number. There is a considerable variation in the evolution of the modal coefficients as the time proceeds. The Galerkin projection is unbounded with less number of modes and gives nonphysical results for such complex flows. In order to recover the correct physics using Galerkin projection, we will have to use an increased number of modes and this will lead to increased computational cost. However, we are interested in recovering the accurate physics as much as possible with less computational cost.

Figure 15 presents results for all DNN frameworks with $p=1$ in the input training data for Ra = $9.4\times 10^{5}$. We can see that the DNN prediction is bounded for all DNN frameworks. Although there exist necessary and sufficient conditions for global boundednessSchlegel and Noack (2015) for Galerkin systems, we have not performed rigorous analysis about boundedness of the proposed frameworks. A deeper discussion of a bounded-input-bounded-output systems can be found elsewhereShi and Han (2007). We observe that the DNN-R and DNN-B frameworks perform better than DNN-S framework especially in the in-sample zone (i.e., $t=0$ to $t=50$). However, we see that the DNN-R and DNN-B frameworks overpredict both modal coefficient $a_{1}$ and $a_{7}$ for out-of-sample zone than the true projection of FOM solution on reduced order space. Figure 16 shows the similar results for all DNN frameworks with $p=4$ in the input training data. We notice that the prediction has improved for all DNN frameworks when we use the modal coefficient history in the input training data. All DNN frameworks are almost able to predict the true projection of FOM solution accurately in the in-sample zone with $p=4$. The problem of overprediction is also reduced for DNN-R and DNN-B frameworks by including the time history in the input data.

Interestingly, including time history in the input training data seems to have helped neural networks to predict more accurate results for both cases (Ra = $3.4\times 10^{5}$ and Ra = $9.4\times 10^{5}$). One explanation for this behavior is that short-term past history of the system helps the neural network to learn the state of the system and predict future state more accurately. Our numerical experiments show that increasing the past history of the system might not help beyond some point. In some cases, including long-term history might give adverse results due to overfitting, or DNN trying to find an unrelated pattern between the input and the output.

Figures 13-16 show the vorticity modal coefficient only for two modes $a_{1}$ and $a_{7}$. We use total RMSE given by Equation (6) to measure the quantitative performance for all DNN frameworks. The total RMSE measures the deviation between the true projection of FOM solution and the prediction by the neural network for all modes. In Table 3, we report the total RMSE for all DNN frameworks investigated in this study for vorticity and temperature modal coefficients for both Rayleigh numbers Ra = $3.4\times 10^{5}$ and Ra = $9.4\times 10^{5}$. Table 3 also reports the root mean square error for ROM-G framework. It can be easily seen that all DNN frameworks perform better than ROM-G framework for both cases. At higher Rayleigh number, the ROM-G framework is unbounded and hence the error is very large. For both cases, we see an improvement in prediction in terms of RMSE for all DNN frameworks as we increase $p$ from $1$ to $4$.

After comparing the time evolution of vorticity, and temperature modal coefficient for all DNN frameworks, we proceed to compare the performance of proposed DNN frameworks in predicting the temperature field in the cavity at final time $t=100$. We compare our results with the true projection of FOM solution on the reduced order space. The FOM solution is obtained by DNS and is discussed briefly in Section IV.2. For lower Rayleigh number case, first 10 modes capture more than $99\%$ of the total energy. Hence, the FOM solution and its true projection will be close to each other. However, for higher Rayleigh number case, first 10 modes capture only 69% of total energy and we expect to see some discrepancy between FOM solution and its true projection on reduced order space. We train the neural network using the evolution of modal coefficients for the true projection of FOM solution and hence we can recover at most true projection of FOM solution.

We compute the instantaneous temperature field using Equation (38). The coefficient $b_{k}$ is considered at the final time step $t=100$. Figure 17 displays the temperature field at the final time predicted by FOM simulation, true projection of FOM on reduced order space, ROM-G framework, and all DNN frameworks. We observe that the FOM solution and its true projection are almost identical. We see that there is some deviation in the temperature field predicted using ROM-G framework and the true projection. The temperature field predicted using DNN-S and DNN-B frameworks are very close to the true projection solution. We see some variation in the solution predicted by DNN-R framework, which can be attributed to the phase shift in modal coefficient predicted by DNN-R framework with $p=1$ (as can be seen in Figure 13). Figure 18 shows the similar results for lower Rayleigh number case with $p=4$. We notice an improvement in prediction by the DNN-R framework and the predicted temperature field is close to the true projection solution.

Figure 19 illustrates the similar results for temperature field at final time $t=100$ for higher Rayleigh number case. The neural network is trained using $p=1$ in the input training data. We see some of the discrepancies between the FOM solution and its projection on reduced order space. The discrepancy is due to less amount of energy in the first 10 POD modes, and hence, some of the flow features get neglected. The deviation is mainly observed at bottom and top of the heated cavity due to the vortical structures formed in these regions at higher Rayleigh number. From Figure 19 we see that the solution predicted by ROM-G framework is very different from the true projection for higher Rayleigh number case. The solution predicted by all DNN frameworks is not identical to the true projection solution. The difference is primarily seen at the bottom and top region. However, we see overall good qualitative agreement between the temperature field predicted by DNN and true projection solution. Figure 20 shows similar results for higher Rayleigh number case when the neural network is trained using $p=4$ in the input training data. We observe a slight improvement in the temperature field prediction. This is consistent with an improvement in the prediction of modal coefficient (refer to Figure 16) with an increase in the solution history in the input data for the neural network. The prediction of modal coefficients is very close to the modal coefficients for the true projection of FOM solution in the in-sample zone ($t=0$ to $t=50$). We do not get the similar accuracy for the out-of-sample zone ($t=50$ to $t=100$). Therefore, we are not able to recover the true projection solution exactly. Despite its limitations, it is clear that our data-driven non-intrusive framework is robust and accurate compared to the intrusive ROM-G framework, and can be used for challenging flow problems where the flow is turbulent and not uniform.

Next, we calculate the time-averaged Nusselt number along the left wall ($x=0$). This quantity can be of interest in many engineering applications. The instantaneous Nusselt number is calculated using Equation (25). The temperature in Equation (25) can be obtained using Equation (38). The gradient of the temperature is computed using a right-sided second-order finite difference scheme. The integration of the temperature gradient is computed using Simpson’s integration rule. After calculating the instantaneous temperature at every time step, we take its average to get the time-averaged Nusselt number. The detail formulas for calculating the instantaneous Nusselt number are provided in the appendix.

In Table 4, we list the statistics of time-averaged Nusselt number for FOM solution, true projection of FOM solution, ROM-G framework, and all DNN frameworks investigated in this study. We can notice all DNN frameworks and intrusive ROM-G framework gives an accurate prediction of the time-averaged Nusselt number. The standard deviation of Nusselt number is also predicted correctly by all DNN frameworks and ROM-G framework for lower Rayleigh number case. We do not recover similar results for higher Rayleigh number case. The time-averaged Nusselt number for true projection solution is slightly different from the FOM solution. The ROM-G framework overpredicts the time-averaged Nusselt number and the difference between the true projection solution and ROM-G prediction is significant. On the other hand, all DNN frameworks predicted the time-averaged Nusselt number close to the true projection solution. The standard deviation of the Nusselt number is also predicted with sufficient accuracy for all DNN frameworks.

To illustrate the temporal variation of Nusselt number, we show the evolution of Nusselt number for FOM, true projection, ROM-G, and DNN-R framework for both Rayleigh numbers in Figure 21. It can be clearly seen that the ROM-G framework fails to predict the temporal behavior of the Nusselt number correctly especially at Ra = $9.4\times 10^{5}$. At lower Rayleigh number, Ra = $3.4\times 10^{5}$, there is a phase shift in the Nusselt number prediction by the ROM-G framework. At this Rayleigh number, with limit cycle oscillations, the Nusselt number is slightly underpredicted (up to second digit accurate) by the true projection of FOM solution on reduced order bases. All DNN frameworks predict the Nusselt number accurately close to the true projection results. We would like to again emphasize that the neural network is trained using the true projection of FOM solution and hence we can at most recover the true projection results and not the FOM results. We also note that much simpler methods, like autocorrelation analysis (instead of the heavy DNN), can be used in predicting this stationary time series problem Shumway and Stoffer (2017). For higher Rayleigh number, the Nusselt number prediction for the ROM-G framework becomes unbounded in the beginning and then it calculates the overpredicted value of instantaneous Nusselt number similar to modal coefficients for vorticity and temperature field. The DNN-R framework correctly predicts the instantaneous Nusselt number close to the true projection of the FOM solution in the in-sample zone and sufficiently accurate results for the out-of-sample zone. To avoid redundancy, we do not present results for other DNN-frameworks since we get a similar prediction.

To summarize, we demonstrated the capability of our DNN frameworks within the non-intrusive ROM setup for the differentially heated cavity problem at two Rayleigh numbers. We do a systematic analysis of our DNN frameworks in terms of prediction of the time evolution of modal coefficients, instantaneous temperature field prediction at the final time, and prediction of time-averaged quantities. Our non-intrusive ROM framework gives sufficiently accurate results for simple as well as complex flows.

VI Concluding Remarks

In this work, we put forward a non-intrusive reduced order modeling framework which uses a deep neural network to predict the dynamics of ROM for complex flow problems. Deep neural networks are capable of approximating a complex nonlinear relationship between the input and output and we achieve this via supervised learning task. The key difference between proposed DNN frameworks is the output variable that is learned by the neural network. The non-intrusive ROM is devised using an encoder-decoder approach and DNN is used for predicting the modal coefficients in an iterative fashion. We use our DNN frameworks with multiple temporal legs (short term history of the state of the system) in the input data. This enables the neural network to take the memory effect into account for predicting the future state of the system. We leverage the classical numerical schemes (backward-difference) in our proposed DNN-B framework, and this framework can also be implemented with other numerical schemes such as Adams-Bashforth, and Adams-Moulton families. First, we use all DNN frameworks for two benchmark problems: Kraichnan-Orszag system and Lorenz system. All DNN frameworks are able to correctly predict each state of the Kraichnan-Orszag system. Even though all DNN frameworks fail to predict correct trajectories for each state of the Lorenz system for a longer duration of time, it predicts the correct dynamics of Lorenz system in terms of the Lorenz attractor.

After evaluating all DNN frameworks for predicting the dynamics of nonlinear dynamical systems, we proceed to reduce order modeling of the differentially heated cavity problem. We extract 1000 snapshots from DNS simulation after the steady state has been reached. The POD bases are constructed using these 1000 snapshots. Based on existing literature and our findings from DNS simulation, we see the change from periodic flow to turbulent flow with an increase in Rayleigh number. For this reason, we test our proposed non-intrusive ROM framework for two cases: lower Rayleigh number case (Ra = $3.4\times 10^{5}$) and higher Rayleigh number case (Ra = $9.4\times 10^{5}$). We use 10 POD modes for our analysis. Due to the turbulent nature of flow at higher Rayleigh number, the first 10 modes capture only 69% of the energy. We sacrifice the advantage of ROM if we increase the number of modes ($R=40$ for 95% of the energy) and hence we attempt to get results close to the true projection of FOM solution on reduced order space with the proposed non-intrusive ROM framework. We assess the performance of non-intrusive ROM framework using different parameters such as prediction capability of modal coefficients, prediction of instantaneous temperature at the final time, and prediction of engineering quantities of interest such as time-averaged Nusselt number. The non-intrusive ROM frameworks perform exceptionally well for all these parameters for lower Rayleigh number case. We see some deviation with the prediction of modal coefficients for higher Rayleigh number case (especially for the out-of-sample zone). Despite this deviation, the proposed framework is able to predict the instantaneous temperature field and time-averaged Nusselt number with sufficient accuracy.

We also compared our results with the results for Galerkin projection (ROM-G framework). The ROM-G framework gives a good prediction for low Rayleigh number case. However, the ROM-G framework is unbounded for higher Rayleigh number case, and produces the wrong prediction. Our analysis for differentially heated cavity problem at two Rayleigh numbers indicates that the non-intrusive ROM setup equiped with the DNN frameworks yields satisfactorily accurate results, and has a potential for reduced order modeling of complex flow problems. In this paper, we used snapshot POD to represent the high-dimensional data onto low-dimensional space. At this point, we can ask the following question: is the POD preconditioning for the model identification step necessary? Given the analogy between POD and shallow autoencoder, it is possible to use deep neural networks to provide a more compact representation of high-dimensional dataBrunton, Noack, and Koumoutsakos (2019). We might arguably assume that a DNN could also learn the potentially much lower dimensional manifold from snaphsot dataLoiseau, Noack, and Brunton (2018). Indeed, this is an exciting future direction and the work is in progress on this topic.

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. O.S. gratefully acknowledges their support. Direct numerical simulations for this project were performed using resources of the Oklahoma State University High Performance Computing Center. We also gratefully acknowledge the support of NVIDIA Corporation with the donation of the GeForce Titan Xp GPU for our research. 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.

Appendix: Nusselt number calculation

The instantaneous Nusselt number of the flow along the left wall ($x=0$) is given by

[TABLE]

where $H$ is the height of the cavity.

For the FOM simulation, we have the temperature data available for each discrete point at each time step. The gradient of the temperature is computed using the right-sided finite difference scheme along the left wall ($x=0$) as given by the below equation

[TABLE]

where $j$ represents the discrete spatial location in $y$-direction, and $h$ is the grid spacing in $x$-direction. The numerical integration is performed using the fourth-order accurate Simpson’s rule. The mean and the standard deviation of the Nusselt number are then computed from the series of instantaneous Nusselt number evaluated at each time step between $t=0$ to $t=100$ using below formulas

[TABLE]

where $N$ is the total number of time steps between $t=0$ to $t=100$.

In the ROM framework, the full order solution can be recovered using the modal coefficient using Equation (38) and is also given below

[TABLE]

where $R$ is the number of modes retained for POD ($R=10$ in this study).

We can compute the temperature field at each time step using Equation (A5) and then follow the same procedure similar to FOM solution. Another faster method to determine the Nusselt number is to calculate the gradient of the mean temperature and each basis function and store it in the memory. The stored gradient can be used to find the instantaneous temperature gradient using equation given below

[TABLE]

If we substitute Equation (A6) in Equation (A1) we get

[TABLE]

where the predetermined coefficients are

[TABLE]

The gradient for the mean temperature and basis function is evaluated using Equation (A2). The temporal value of temperature modal coefficient is determined using Galerkin projection for ROM-G framework and using DNN frameworks for non-intrusive ROM setup. The numerical integration of instantaneous temperature gradient along the left wall is computed using the composite Simpson’s rule with SciPy function integrate.simps for all ROMs. The mean and standard deviation of the Nusselt number is then computed using equations (A3) and (A4), respectively.