Closure Learning for Nonlinear Model Reduction Using Deep Residual Neural Network

Abstract

Developing accurate, efficient, and robust closure models is essential in the construction of reduced order models (ROMs) for realistic nonlinear systems, which generally require drastic ROM mode truncations. We propose a deep residual neural network (ResNet) closure learning framework for ROMs of nonlinear systems. The novel ResNet-ROM framework consists of two steps: (i) In the first step, we use ROM projection to filter the given nonlinear system and construct a spatially filtered ROM. This filtered ROM is low-dimensional, but is not closed. (ii) In the second step, we use ResNet to close the filtered ROM, i.e., to model the interaction between the resolved and unresolved ROM modes. We emphasize that in the new ResNet-ROM framework, data is used only to complement classical physical modeling (i.e., only in the closure modeling component), not to completely replace it. We also note that the new ResNet-ROM is built on general ideas of spatial filtering and deep learning and is independent of (restrictive) phenomenological arguments, e.g., of eddy viscosity type. The numerical experiments for the 1D Burgers equation show that the ResNet-ROM is significantly more accurate than the standard projection ROM. The new ResNet-ROM is also more accurate and significantly more efficient than other modern ROM closure models.

1. Introduction

Many scientific and engineering applications, such as weather forecasting, ocean modeling, and cardiovascuar flow simulation, can often be represented by multiscale systems of ordinary differential equations (ODE) in high-dimensional modal spaces. The analysis and high-fidelity simulation of such systems can be very expensive even on high-performance computing systems. Consequently, using the full order model (FOM) for such simulations can be impractical for time critical applications, such as flow control and parameter estimation. To alleviate the computation burden of the FOM simulation, reduced order models (ROMs) have been successfully used.

ROMs seek a low-dimensional approximation of the FOM with orders of magnitude reduction in computational cost. The classical projection based ROM approach is first to construct a low dimensional space using data-driven reduction methods, such as proper orthogonal decomposition (POD) or dynamical modal decomposition (DMD). The ROM dynamics can be obtained via Galerkin projection of the FOM onto the reduced space [1,2,3,4,5].

The Galerkin projection reduced order model (GP-ROM) can be very efficient and relatively accurate for many systems [3,4]. There are, however, systems (e.g., convection-dominated fluid flows) for which the GP-ROM can generate inaccurate approximations. There are several approaches to address this inaccuracy (see, e.g., [6,7,8,9]). In this paper, we focus on one of the main reasons for this inaccuracy: the ROM closure problem, i.e., modeling the interaction between the GP-ROM modes and the discarded modes. Indeed, due to the inherently drastic mode truncation required in realistic settings, the dimension of the GP-ROM space is too low to resolve the complex nonlinear interactions of the fluid system [1,2,10,11]. GP-ROMs that do not include a ROM closure model can yield inaccurate results, often in the form of spurious numerical oscillations [1,2,12,13]. Endowing GP-ROM with closure models could extend the applicability of GP-ROM in many fluid mechanics applications, such as flow control, climate modeling, and weather forecasting [1,2].

ROM closure models for nonlinear systems have been proposed in, e.g., [14,15,16,17,18,19,20,21,22,23,24]. The vast majority of the current ROM closure models aim at mitigating the numerical instability observed in GP-ROMs that do not include a closure model. Some of these ROM closure models use stabilization techniques that have been developed in standard discretization methods (e.g., in the finite element community) [25,26,27]. Other ROM closure models have imported ideas developed in standard CFD methodologies, e.g., large eddy simulation (LES) [1,24,28,29]. The overwhelming majority of the current ROM closure models can be categorized as stabilization techniques (for a notable exception, see the approximate deconvolution ROM closure model [28] that uses a mathematical framework inspired from image processing).

This is in stark contrast with classical LES, where a plethora of closure models have been proposed over the years. The main difference between ROM closure and LES closure is that the latter has been entirely built around physical insight from the statistical theory of turbulence (e.g., energy cascade and Kolmogorov’s theory), which is generally posed in the Fourier setting [30,31]. This physical insight is generally not available in the ROM setting (see, however, [32] for initial steps). Thus, current ROM closure models have generally been deprived of this powerful methodology that represents the core of most LES closure models.

We believe that machine learning represents a natural solution for ROM closure modeling. Indeed, since physical insight cannot be easily extended to the ROM setting, available data and machine learning can be utilized instead to develop ROM closure models.

We propose a novel ROM closure modeling framework that is constructed by using available data and deep residual neural network (ResNet) (for details, see, e.g., [33,34,35,36,37]). The resulting ROM, which we call the residual neural network ROM (ResNet-ROM), is schematically illustrated in (1) and Figure 1 (see Section 3 for details). We emphasize that, in the new ResNet-ROM framework, data is used only to complement classical physical modeling (i.e., only in the closure modeling component) [29,38], not to completely replace it [39]. Thus, the resulting ResNet-ROM combines the strengths of both physical and data-driven modeling.

$$\begin{array}{ccccc}\mathrm{FOM}& \stackrel{\mathrm{filtering}}{\to }& \mathrm{GP-ROM}& \stackrel{\mathrm{learning\; closure}}{\to }& \mathrm{ROM}\end{array}$$

(1)

The main contributions of this paper can be summarized as follows:

• A novel ROM closure learning framework centered around deep neural networks.
• A hybrid framework that synthesizes the strengths of physical modeling and data-driven modeling.
• Very good performance in numerical tests, in both the reconstructive and the predictive regime.
• Significant improvement in numerical accuracy compared with state of the art ROM closure models.

Machine learning has recently been utilized to construct ROM closure models (see, e.g., [22,38,40,41,42]). The ROM closure model proposed in [22] is similar to the new ResNet-ROM. There are, however, two major differences between the two ROM closure models: The first difference is that the ROM closure model in [22] uses FOM data to model both the linear and the nonlinear terms in the underlying equations (see Section 5.1 in [22]). In contrast, in the LES spirit [31], the ResNet-ROM models only the nonlinear terms (see Algorithm 1). One motivation for modeling only the nonlinear terms is given in the theoretical and numerical investigation in [43], where it was shown that the contribution of the linear terms to the ROM closure term (i.e., the commutation error) is not significant in convection-dominated problems, such as those we consider in this study. The second major difference between the ROM closure model in [22] and the ResNet-ROM is the machine learning approach used to construct the ROM closure model: Bayesian regularization and extreme learning machine approaches are used in [22], whereas a deep residual neural network (ResNet) is used to construct the ResNet-ROM.

2. Reduced Order Model

To construct the novel ResNet-ROM framework, we use the Navier-Stokes equations (NSE) for incompressible fluid flows:

$$\left\{\begin{array}{cc}\hfill \frac{\partial \mathit{u}}{\partial t}-R{e}^{-1}\mathsf{\Delta }\mathit{u}+(\mathit{u}·\nabla )\mathit{u}+\nabla p& =0,\hfill \\ \hfill \nabla ·\mathit{u}& =0,\hfill \end{array}\right.$$

(2)

which are defined on a spatial domain $\mathsf{\Omega }$ and on the time interval $[0,T]$. In the NSE (2), $\mathit{u}$ is the velocity, p the pressure, and $Re$ the Reynolds number. The ROM basis $\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_r\}$, where r is small, represents the large, energy containing structures in the flow, and is obtained by using available numerical or experimental data and, e.g., the POD or DMD methods [1,3,4]. The ROM velocity approximation is defined as

$${\mathit{u}}_r(\mathbf{x},t)\equiv \sum _{j=1}^r{a}_j\left(t\right){\mathit{\phi }}_j\left(\mathbf{x}\right),$$

(3)

where ${\left\{a_j\left(t\right)\right\}}_{j=1}^r$ are the sought time-varying coefficients. To determine these coefficients, we use a Galerkin procedure. We replace $\mathit{u}$ in (2) with the ROM approximation (3), we take the $L^2$ inner product between the resulting system and each ROM basis function $\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_r\}$, and then we integrate by parts: $\forall i=1,\dots ,r,$

$$\begin{array}{}\hfill \left(\frac{\partial {\mathit{u}}_r}{\partial t},{\mathit{\phi }}_i\right)+\frac{1}{Re}\left(\nabla {\mathit{u}}_r,\nabla {\mathit{\phi }}_i\right)+\left(({\mathit{u}}_r·\nabla ){\mathit{u}}_r,{\mathit{\phi }}_i\right)=0.\end{array}$$

To derive the above equation, we assumed that the ROM velocity modes are perpendicular to the discrete pressure space, which is the case if standard mixed FEs (e.g., Taylor-Hood) are used for the snapshot creation [1,29]. Using (3), we obtain the standard Galerkin projection ROM (GP-ROM):

$$\begin{array}{c}\hfill \dot{\mathit{a}}=A\mathit{a}+{\mathit{a}}^{\top }B\mathit{a},\end{array}$$

(4)

which can be written componentwise as follows: $\forall i=1\dots r,$

$${\dot{a}}_i=\sum _{m=1}^r{A}_{im}{a}_m\left(t\right)+\sum _{m=1}^r\sum _{n=1}^r{B}_{imn}{a}_n\left(t\right)a_m\left(t\right),$$

(5)

where $A_{im}=-R{e}^{-1}\left(\nabla {\mathit{\phi }}_m,\nabla {\mathit{\phi }}_i\right)$ and $B_{imn}=-\left({\mathit{\phi }}_m·\nabla {\mathit{\phi }}_n,{\mathit{\phi }}_i\right)$.

3. Closure Learning

3.1. Residual Neural Network (ResNet)

Deep residual neural network (ResNet) has been first introduced for image recognition in [33]. ResNet has been widely studied and applied in many supervised learning tasks. Recent mathematical understanding of deep ResNet has been achieved in the ODE representation of ResNet; for a comprehensive introduction, see, e.g., [35,36,37,44].

To construct the novel ResNet-ROM framework, we consider the ResNet model, which is illustrated in Figure 2.

The input values of forward propagation in the ResNet are given by

$$X_{t+1}=X_t+tanh(W_t{X}_t+b_t),t=1,\dots ,N-1,$$

(6)

where N is the number of layers in the network architecture, and $X_t\in {\mathbb{R}}^s$ is the output from each ResNet block at time t. $W_t$ and $b_t$ are the weight matrix and bias at each layer, respectively. The ResNet propagation starts from step $t=1$ with the nonlinear activation function tanh. The initial input layer for the network is $X_0={[Re,t,a_1,\dots ,a_r]}^T$. For a standard ResNet in image classification, a convolution layer (CNN) is often included in a residual block. In our work, we use a simplified version of the ResNet that does not include CNN.

3.2. ROM Closure Modeling

In realistic nonlinear systems (e.g., convection-dominated flows), current GP-ROMs of the form (4) generally yield inaccurate results, often in the form on numerical instabilities. This inaccurate behavior is due to the fact that, in order to maintain a low computational cost, GP-ROMs are constructed with a drastically truncated ROM basis $\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_r\}$, where r is generally small (e.g., $r=\mathcal{O}\left(10\right)$). In realistic nonlinear systems, this extreme truncation cannot capture the complex nonlinear interactions among the various degrees of freedom of the system. Thus, current GP-ROMs are computationally efficient, but numerically inaccurate. To alleviate this inaccurate behavior, GP-ROMs are often supplemented with a ROM closure model [1,24], i.e., a model for the important interactions between the ROM basis $\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_r\}$ and the discarded ROM modes $\{{\mathit{\phi }}_{r+1},\dots ,{\mathit{\phi }}_R\}$, where R is the dimension of the input data. For example, for the NSE, the standard GP-ROM (4) is generally modified as follows:

$$\begin{array}{c}\hfill \dot{\mathit{a}}=A\mathit{a}+{\mathit{a}}^{\top }B\mathit{a}+\mathsf{\tau },\end{array}$$

(7)

where $\mathsf{\tau }$ is a ROM closure model that represents the interactions between the ROM modes and the discarded modes.

We emphasize that the same closure problem needs to be addressed when classical numerical discretization schemes (e.g., finite element or spectral methods) are used in the numerical simulation of turbulent flows. In those settings, classical discretizations schemes are used in inherently under-resolved regimes (i.e., they use coarse meshes or too few spectral modes). To ensure the relative accuracy of these under-resolved simulations, various types of closure models for the unresolved (e.g., subgrid-scale) information are generally utilized. These closure models are central to, e.g., large eddy simulation (LES) [31], one of the main approaches to the numerical simulation of turbulent flows. The vast majority of LES closure models have been constructed by using physical insight from Kolmogorov’s statistical theory of turbulence. The concept of energy cascade is central in the development of LES closure models. The energy cascade states that energy enters the system at the large scales, is transferred to smaller and smaller scales through nonlinear interactions, and is dissipated at the smallest scale (i.e., the Kolmogorov’s scale). Thus, most LES closure models (e.g., of eddy viscosity type) aim at recovering the energy cascade displayed by the original system (i.e., the NSE).

This physical insight cannot be easily extended to the ROM setting (see, however, [32] for a preliminary numerical investigation). Thus, current ROM closure models have generally been deprived of this powerful methodology that represents the main tool in the development of most LES closure models.

3.3. ROM Closure Learning

Our vision is that machine learning represents a natural approach for constructing ROM closure modeling. Indeed, since physical insight cannot generally be used in a ROM setting, data and machine learning can be utilized instead to develop ROM closure models. Furthermore, data is used only to construct the ROM closure model; the other ROM operators are built by using the classical Galerkin projection. Thus, data and machine learning complement (instead of replace) physical based modeling, yielding a hybrid framework that synthesizes the strengths of both approaches [17,18,22,27,29,42,45,46,47,48].

In this section, we propose a novel ROM closure modeling framework that is constructed by using available data and the ResNet approach described in Section 3.1. The resulting ROM, which we call the residual neural network ROM (ResNet-ROM), is schematically illustrated in (1) and Figure 1, and is summarized in Algorithm 1.

To obtain the explicit formula (11) in Algorithm 1, we develop a large eddy simulation ROM (LES-ROM) framework [24,28]: First, we filter the high-resolution (R-dimensional) ROM approximation of the NSE with a low-pass ROM spatial filter, denoted by overbar in (11). In this paper, as a ROM spatial filter, we use the ROM projection onto the space spanned by the first r ROM modes [29,48] (see [28] for an alternative ROM spatial filter): For a given $u\in \mathrm{span}\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_R\}$, the ROM projection seeks $\overline{u}\in \mathrm{span}\{{\mathit{\phi }}_1,\dots ,{\mathit{\phi }}_r\}$ such that

$$\begin{array}{c}\hfill (\overline{u},{\phi }_i)=(u,{\phi }_i)\forall i=1,\dots ,r,\end{array}$$

(8)

where $(·,·)$ denotes the $L^2$ inner product. The resulting filtered equations can be cast within a variational multiscale (VMS) framework [49], yielding a two-scale VMS-ROM [50]. Equation (11) represents the explicit VMS-ROM formula for the ROM closure term $\mathsf{\tau }$ in(10). In the two-scale VMS-ROM [50], the ROM closure model is developed in two steps. First, an ansatz is used to approximate the ROM closure term: $\mathsf{\tau }\approx \tilde{A}\mathit{a}+{\mathit{a}}^{\top }\tilde{B}\mathit{a}$. Then, FOM data for $\mathsf{\tau }$ in the training time interval is used to solve a least squares problem to determine the optimal entries in $\tilde{A}$ and $\tilde{B}$. In the new ResNet-ROM, we construct the ROM closure model in a fundamentally different way. Instead of making an ansatz on the structural form of the ROM closure term $\mathsf{\tau }$ (as in the two-scale VMS-ROM [50]), in Algorithm 1 we consider a general structural formulation and use the ResNet approach to construct the ROM closure model. To this end, in (6), we make the following choices: The initial layer contains the Reynolds number ($Re$), the current time (t), and the current ROM coefficients in (3) ($a_1,\dots ,a_r$), i.e., $X_0=[Re,t,a_1,\dots ,a_r]\in {\mathbb{R}}^{(s+2)}$. The final output layer provides an approximation for the closure term, i.e., $X_N\approx {\mathsf{\tau }}^{FOM}$. The optimization problem associated with this ResNet is given by

$$min\parallel {\mathsf{\tau }}^{ansatz}-{\mathsf{\tau }}^{FOM}{\parallel }_F^2+\lambda R(W,b),$$

(9)

where the regularizer R penalizes undesirable parameters and can prevent overfitting, ${\mathsf{\tau }}^{ansatz}=f(\mathit{a},Re,t_j)$ is the output from the neural network, W and b are weights in the network, ${\parallel ·\parallel }_F$ is the Frobenius norm, and $\lambda$ is a hyperparameter for the $L^2$ regularization to prevent overfitting [36,37]. The minimization problem (9) is over the parameters defining the structural form of the function f used to define ${\mathsf{\tau }}^{ansatz}$.

Algorithm 1 ResNet-ROM

1: Consider the ROM closure model $$\begin{array}{c}\hfill \dot{\mathit{a}}=A\mathit{a}+{\mathit{a}}^{\top }B\mathit{a}+\mathsf{\tau }.\end{array}$$ (10) 2: Use snapshot data to compute the true vector $\mathsf{\tau }$ in (10), ${\mathsf{\tau }}^{FOM}$: $$\begin{gathered} \begin{array}{}\hfill {\tau }_i^{FOM}\left(t_j\right)& =& -({\overline{\left({\mathit{u}}_R^{FOM}\left(t_j\right)·\nabla \right){\mathit{u}}_R^{FOM}\left(t_j\right)}}^r\hfill \\ \\ -\left({\mathit{u}}_r^{FOM}\left(t_j\right)·\nabla \right){\mathit{u}}_r^{FOM}\left(t_j\right),{\mathit{\phi }}_i),\end{array} \end{gathered}$$ (11) where an overbar indicates spatial filtering with ROM projection. 3: Use snapshot data and (11) to define the approximation function $\mathsf{\tau }$ in (10), ${\mathsf{\tau }}^{ansatz}$: $$\begin{array}{c}\hfill {\mathsf{\tau }}^{ansatz}\left(t_j\right)=f(\mathit{a},Re,t_j),\end{array}$$ (12) where f is a generic function that needs to be determined. 4: Use ResNet to train the closure term, i.e., to find the form of f in(12) that is optimal with respect to the minimization problem (9). 5: The novel ResNet-ROM has the following form: $$\dot{\mathit{a}}=A\mathit{a}+{\mathit{a}}^{\top }B\mathit{a}+f(\mathit{a},Re,t).$$ (13)

4. Numerical Experiments

4.1. Implementation

As a test problem for our new approach, we use the 1D Burgers equation, which has been used to test new ROM ideas in simplified settings [29,51,52]:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{u}_t-R{e}^{-1}{u}_{xx}+u{u}_x=0\hfill & x\in \mathsf{\Omega },\hfill \\ u(x,0)=u_0\left(x\right)\hfill & x\in \mathsf{\Omega },\hfill \\ u(x,t)=0\hfill & x\in \partial \mathsf{\Omega },\hfill \end{array}\right.\end{array}$$

(14)

where, for consistency with the notation used for the NSE, the diffusion parameter is denoted as $R{e}^{-1}$. In our numerical tests, $\mathsf{\Omega }=[0,1]$ is the computational domain and the time domain is $[0,1]$. We use the same initial conditions as those utilized in [29,51,52]: $u_0\left(x\right)=1,x\in (0,1/2]$, $u_0\left(x\right)=0,x\in (1/2,1]$. These initial conditions yield a steep internal layer that is difficult to capture in the convection-dominated regime that we consider [29,51,52]. We first use a piecewise linear finite element discretization to generate the FOM solution. To this end, we utilize a uniform mesh with $N=1024$ grid points (which yields a meshsize $h=1/1024$) and the forward Euler method with $\mathsf{\Delta }t={10}^{-4}$ for the time discretization. To construct the ROM basis, we collect 101 snapshots sampled from $[0,1]$. We build the neural network in PyTorch and we train the ROM closure model with a 6 block ResNet with the common Adam optimizer [53]. We perform all the computational experiments on a Linux system laptop with Nvidia Geforce GTX GPU hardware.

4.2. Reconstruction

In this section, we consider the reconstructive regime, i.e., we test the ROMs at the same $Re$ as the $Re$ at which the ROMs are constructed. We choose $Re=1000$ in (14) and we use $r=6$ basis functions in all ROMs. In Figure 3, we plot the solutions of the FOM (top left), GP-ROM (top right), and ResNet-ROM (bottom left). When compared with the FOM data, the ResNet-ROM solution is significantly more accurate than the standard GP-ROM solution. In Figure 3, we also plot the FOM, GP-ROM, and ResNet-ROM solutions at the final time step (i.e., at $t=1$). This plot shows that the closure term in the ResNet-ROM plays an important role in stabilizing the ROM approximation. Indeed, the GP-ROM solution displays large, spurious numerical oscillations. These oscillations are dramatically decreased in the ResNet-ROM solution.

In Figure 4, we plot the the time evolution of the ROM coefficients a for the FOM, GP-ROM, and ResNet-ROM. The plots show that the ResNet-ROM is significantly more accurate than the standard GP-ROM.

4.3. Prediction

To study the robustness of the new ResNet-ROM, we test its predictability, i.e., we train the ResNet-ROM closure term on data from multiple $Re$ and we then test the ResNet-ROM to predict the ROM dynamics at different $Re$. The training data space is sampled at $Re=[20,50,100,200,500,800,1000]$, and the test data contains $Re=[30,80,300,1200]$. Note that the solution of the Burgers equation is affected by $Re$. Small $Re$ values yield a slow movement of the sharp internal layer, while large $Re$ values can speed up this movement.

In Figure 5, for $Re=30,80$, and 1200 (which are different from the training $Re$ values), we plot the solutions for the DNS (first column), GP-ROM (second column), ResNet-ROM (third column), and final time solution for all three simulations (fourth column). These plots show that the ResNet-ROM is consistently the most accurate ROM, especially for the largest $Re$ value. In Figure 5, we also plot the FOM, GP-ROM, and ResNet-ROM solutions at the final time step (i.e., at $t=1$). These plots show that the closure term in the ResNet-ROM plays an important role in stabilizing the ROM approximation.

In Figure 6, we plot the the time evolution of the ROM coefficients a for the FOM, GP-ROM, and ResNet-ROM. These plots show that the ResNet-ROM is significantly more accurate than the standard GP-ROM.

Overall, we draw the same conclusion as in the reconstructive regime (Section 4.2): In all cases, the ResNet-ROM is significantly more accurate than the standard GP-ROM in the predictive regime.

4.4. Comparison

In this section, in the numerical simulation of the Burgers equation, we also make a numerical comparison between the new ResNet-ROM and several modern ROM closure models: the POD artificial viscosity model (POD-AV) [51], the POD Smagorinsky model (POD-L) [54], the evolve-then-filter ROM (EF-ROM) [55], the approximate deconvolution ROM (AD-ROM) [28], and the data-driven filtered ROM (DDF-ROM) [29].

In

Table 1. $L^2$ errors of the new ResNet-ROM and other closure models.

Model	${\mathit{L}}^2$ Error
ResNet-ROM	4.32 × 10${}^{-4}$
POD-AV [51]	9.01 × 10${}^{-3}$
POD-L [54]	1.73 × 10${}^{-2}$
EF-ROM [55]	6.99 × 10${}^{-2}$
AD-ROM [28]	6.33 × 10${}^{-2}$
DDF-ROM [29]	6.27 × 10${}^{-2}$

, we list the $L^2$ errors of the ResNet-ROM and the other closure models. These results show that the ResNet-ROM error is at least an order of magnitude lower than the errors of the other ROM closure models.

For the online ResNet-ROM integration, we use the scipy package with the built-in integration function odeint. In the online stage, the new ResNet-ROM is more efficient than the other ROM closure models: The online CPU time of the ResNet-ROM is $3.89$ s, whereas the online CPU times of the EF-ROM, AD-ROM, and DDF-ROM are $6.91,7.26,$ and $4.42$ s, respectively [28,29,55].

The CPU time of the offline training of the new ResNnet-ROM is $122.74$ s for the current dataset with 10,000 epochs (iterations). Thus, even though the online cost for ResNet-ROM is lower than the online cost of the other closure models that are compared in this paper, the neural network training cost of ResNet-ROM is much higher than the offline training cost of the other ROM closure models [28,29,55].

4.5. Sensitivity

In this section, we perform a sensitivity study of the new ResNet-ROM with respect to two parameters: (i) the hyperparameter $\lambda$ used in the regularization of the neural network training; and (ii) the meshsize h used in the snapshot generation. We also investigate the potential improvement in GP-ROM accuracy when the number of snapshots and the dimension r of the ResNet-ROM are increased.

The parameter $\lambda$ is an $L^2$ regularization parameter used in the neural network training to prevent overfitting. In our numerical simulations, we pick the parameter $\lambda$ based on the validation error performance in the training phase. In our training dataset for $Re=[20,50,100,200,500,800,1000]$, we test the following parameter values: $\lambda =0,1,{10}^{-1},{10}^{-2},{10}^{-3}$. In Figure 7, we plot the ResNet-ROM error for different parameter values. The plot in Figure 7 shows that the ResNet-ROM is not very sensitive to the hyperparameter $\lambda$. Thus, in our numerical tests, we fix $\lambda =0.01$ with 10,000 epochs (iterations) in the Resnet-ROM training.

We also perform a sensitivity study for the new ResNet-ROM with respect to the meshsize h. In Figure 8, for the reconstructive regime and $Re=1000$, we plot the FOM, GP-ROM, and ResNet-ROM results for three coarse meshsize values. The plots in Figure 8 show that the ResNet-ROM is still significantly more accurate than the GP-ROM. In the predictive regime, however, both the ResNet-ROM and the GP-ROM yield inaccurate results for these three coarse meshes. The relationship between the FOM mesh utilized to generate the snapshots and the ROM accuracy is subject of current research (see, e.g., [56,57]) and should be further investigated for the new ResNet-ROM.

We also investigate the ResNet-ROM rate of convergence with respect to the meshsize h. To this end, we fixed the ResNet-ROM dimension at $r=30$ and the time-step size at $\mathsf{\Delta }t={10}^{-3}$. The plot in Figure 9 shows that the rate of convergence (obtained with a least squares fit) is about $h^{1.70}$, which is an acceptable approximation to the theoretical rate of convergence of $h^2$.

Finally, we investigate the potential improvement in GP-ROM accuracy when the number of snapshots and the dimension r of the ResNet-ROM are increased.

First, we collect the maximum number of snapshots that are available from the FOM simulations in the training interval. That is, we collect 1001 snapshots, which yield a snapshot matrix whose rank is 251. Comparing the left plot in Figure 10 (for the highest number of snapshots) with the top right plot in Figure 3 (for the lower number of snapshots), we conclude that increasing the number of snapshots does not seem to improve the ResNet-ROM accuracy.

Next, we increase the GP-ROM dimension. In Figure 10, we plot the GP-ROM results for three r values: $r=10$, $r=20$, and $r=30$. As expected, the GP-ROM accuracy improves as we increase r. We emphasize, however, that the ROM closure models (such as that used in the new ResNet-ROM) are designed to improve the GP-ROM accuracy in under-resolved numerical simulations, i.e., when only a few ROM basis functions can be used, which is often the case in realistic settings [1,2,14,16,17,18,20,21,22,23,25,27,28,29,32,38,40,41,42,46,47,48,50,51,54,58]. As shown earlier, for $r=10$, the ResNet-ROM is significantly more accurate than the GP-ROM both in the reconstructive and the predictive regimes.

5. Conclusions

In this paper, we used available data and deep residual neural networks (ResNet) to construct a novel reduced order model (ROM) closure for complex nonlinear settings. We emphasize that the ResNet-ROM closure terms are much more general than the phenomenological ansatzes generally used in ROM closure modeling [29]. We tested the novel ResNet-ROM in the numerical simulation of the Burgers equation. For comparison purposes, we investigated the standard Galerkin projection ROM (GP-ROM) and the full order model (FOM). We considered two settings: (i) a reconstructive regime, in which the Reynolds number $Re$ is the same in the training and testing stages; and (ii) a predictive regime, in which $Re$ used in the testing regime is different from the $Re$ used in the training regime. In both regimes, the new ResNet-ROM was consistently more accurate than the standard GP-ROM. Furthermore, the ResNet-ROM was also dramatically more accurate than several other ROM closure models from the literature.

There are several research directions that we plan to pursue: We will test the novel ResNet-ROM on realistic test problems (e.g., 3D turbulent flows) and we will compare it with state of the art closure models. We will also investigate alternative approaches to develop the ROM closure term $\mathsf{\tau }$. Indeed, in this paper we used the ROM projection as a spatial filter in the construction of the ROM closure term $\mathsf{\tau }$. We plan to investigate different ROM spatial filters, such as the ROM differential filter [55]. This alternative ROM filter will yield a different ROM closure term $\mathsf{\tau }$ and, therefore, a different ResNet-ROM. Finally, a numerical investigation of alternative machine learning approaches (see, e.g., [22]) could yield improved ROM closure models. Indeed, for the one-dimensional Burgers equation test case considered in this paper, the computational cost of training the ROM closure model with ResNet was acceptable. For more complex settings, however, this computational cost could be much higher. In those cases, neural networks with a lower computational cost could be more effective.