Data-Driven Correction Reduced Order Models for the Quasi-Geostrophic Equations: A Numerical Investigation

TL;DR

This study evaluates a data-driven correction reduced order model (DDC-ROM) for simulating quasi-geostrophic equations, demonstrating improved accuracy and stability over traditional methods, especially with constrained variants over extended time intervals.

Contribution

Introduces and tests a data-driven correction ROM with physical constraints, enhancing simulation accuracy and stability for quasi-geostrophic equations.

Findings

DDC-ROM outperforms standard G-ROM in accuracy.

Constrained DDC-ROM further improves stability.

Both models perform well over longer simulation times.

Abstract

This paper investigates the recently introduced data-driven correction reduced order model (DDC-ROM) in the numerical simulation of the quasi-geostrophic equations. The DDC-ROM uses available data to model the correction term that is generally used to represent the missing information in low-dimensional ROMs. Physical constraints are added to the DDC-ROM to create the constrained data-driven correction reduced order model (CDDC-ROM) in order to further improve its accuracy and stability. Finally, the DDC-ROM is tested on time intervals that are longer than the time interval over which it was trained. The numerical investigation shows that, for low-dimensional ROMs, both the DDC-ROM and CDDC-ROM perform better than the standard Galerkin ROM (G-ROM) and the CDDC-ROM provides the best results.

Tables

Table 1. Table 1. The relative errors for the time-averaged streamfunction defined by ( 26 ).

$r$ values	G-ROM	DDC-ROM	CDDC-ROM
$r=5$	n/a	n/a	5.07e+00
$r=10$	2.06e+02	3.25e-01	9.58e-02
$r=15$	3.05e+02	2.83e-01	1.03e-01
$r=20$	5.05e+00	1.39e-01	1.20e-01
$r=25$	1.73e+00	1.61e-01	2.96e-01
$r=30$	1.47e+00	1.18e-01	4.58e-01
$r=35$	6.83e-01	1.34e-01	7.69e-01
$r=40$	4.69e-01	9.17e-02	4.92e-01
$r=45$	2.81e-01	4.16e-02	6.83e+00
$r=50$	3.66e-01	8.20e-02	4.16e-01

Table 2. Table 2. The ROM relative errors for the time-averaged streamfunction defined in ( 26 ) for the two predictive test cases: the POD basis functions are generated using DNS snapshots over the time interval [ 10 , 45 ] 10 45 [10,45] for Case I and over [ 10 , 35 ] 10 35 [10,35] for Case II. The ROM simulations are carried out in the time interval [ 10 , 80 ] 10 80 [10,80] .

$r$ values	Predictive Case I		Predictive Case II
	G-ROM	DDC-ROM	G-ROM	DDC-ROM
$r=5$	n/a	n/a	n/a	n/a
$r=10$	1.39e+04	3.33e-01	1.91e+02	6.56e-01
$r=15$	9.58e+00	3.76e-01	1.94e+01	3.27e-01
$r=20$	5.37e+00	1.35e-01	6.35e+00	1.80e-01
$r=25$	2.28e+00	9.47e-02	1.44e+00	3.05e-01
$r=30$	5.65e-01	1.32e-01	5.27e-01	1.96e-01
$r=35$	2.88e-01	1.76e-01	2.26e-01	1.71e-01
$r=40$	2.07e-01	1.72e-01	2.44e-01	2.20e-01
$r=45$	2.86e-01	2.15e-01	3.20e-01	2.34e-01
$r=50$	1.67e-01	2.09e-01	3.53e-01	2.72e-01

Equations

$$\displaystyle\bigl{(}\overset{\bullet}{\bm{u}}\,,\bm{v}\bigr{)}=\bigl{(}\bm{f}(\bm{u})\,,\bm{v}\bigr{)}\,,\qquad\forall\,\bm{v}\in\bm{X}\,,$$

$$\displaystyle\bigl{(}\overset{\bullet}{\bm{u}_{r}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}\,,\qquad i=1,\ldots,r.$$

$$\displaystyle\bigl{(}\overset{\bullet}{\bm{u}_{R}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\bm{f}(\bm{u}_{R})\,,\bm{\varphi}_{i}\bigr{)}\,,\qquad i=1,\ldots,r.$$

$$\displaystyle\bigl{(}\overset{\bullet}{\bm{u}_{r}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\bm{f}(\bm{u}_{R})\,,\bm{\varphi}_{i}\bigr{)}\,,\qquad i=1,\ldots,r,$$

$$\displaystyle\bigl{(}\overset{\bullet}{\bm{u}_{r}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}+\bigl{[}\bigl{(}\bm{f}(\bm{u}_{R})\,,\bm{\varphi}_{i}\bigr{)}-\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}\bigr{]},\qquad i=1,\ldots,r.$$

$$\displaystyle\boxed{\text{Correction}=\bigl{[}\bigl{(}\bm{f}(\bm{u}_{R})\,,\bm{\varphi}_{i}\bigr{)}-\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}\bigr{]}.}$$

$$\displaystyle\text{Correction}=\bigl{[}\bigl{(}\bm{f}(\bm{u}_{R})\,,\bm{\varphi}_{i}\bigr{)}-\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}\bigr{]}\approx\bigl{(}{\bf g}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}\,,$$

$$\displaystyle\boxed{\bigl{(}\overset{\bullet}{\bm{u}_{r}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\bm{f}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)}+\bigl{(}{\bf g}(\bm{u}_{r})\,,\bm{\varphi}_{i}\bigr{)},}\qquad i=1,\ldots,r.$$

$$\displaystyle\boxed{\min_{{\bf g}\text{ parameters}}\ \sum_{j=1}^{M}\left\|\text{Correction}(t_{j})-\bigl{(}{\bf g}(\bm{u}_{r}(t_{j}))\,,\bm{\varphi}_{i}\bigr{)}\right\|^{2}\,.}$$

$$\frac{\partial \omega }{\partial t}+J(\omega ,\psi )-R{o}^{-1}\frac{\partial \psi }{\partial x}$$

$$\omega$$

$$F=\sin (\pi (y-1)).$$

$$\psi (t,x,y)=0,\omega (t,x,y)=0\text{for}(x,y)\in \partial \Omega \text{ and }t\ge 0.$$

$$\left(f(\omega ),v\right)=-\left(J(\omega ,\psi ),v\right)+R{o}^{-1}\left(\frac{\partial \psi }{\partial x},v\right)-R{e}^{-1}\left(\nabla \omega ,\nabla v\right)+R{o}^{-1}\left(F,v\right).$$

$$X^r:=\text{span}\{{\phi }_1,\dots ,{\phi }_r\},$$

$$-\Delta {\varphi }_i(x,y)={\phi }_i(x,y),\text{ subject to }{\varphi }_i(x,y)=0,\text{for}(x,y)\in \partial \Omega ,i=1,2,\cdots ,r.$$

$$\left(\frac{\partial {\omega }_r}{\partial t},{\phi }_i\right)+\left(J({\omega }_r,{\psi }_r),{\phi }_i\right)-R{o}^{-1}\left(\frac{\partial {\psi }_r}{\partial x},{\phi }_i\right)+R{e}^{-1}\left(\nabla {\omega }_r,\nabla {\phi }_i\right)=R{o}^{-1}\left(F,{\phi }_i\right),$$

$$\dot{\mathbf{a}}=\mathbf{b}+A\mathbf{a}+{\mathbf{a}}^{\top }B\mathbf{a},$$

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

$$=$$

$$\displaystyle\qquad\min_{\tilde{A}}\ \sum_{j=1}^{M}\,\biggl{\|}\biggl{[}\left(\frac{\partial\omega_{R}}{\partial x}\frac{\partial\psi_{R}}{\partial y}-\frac{\partial\psi_{R}}{\partial x}\frac{\partial\omega_{R}}{\partial y},\varphi_{i}\right)-\left(\frac{\partial\omega_{r}}{\partial x}\frac{\partial\psi_{r}}{\partial y}-\frac{\partial\psi_{r}}{\partial x}\frac{\partial\omega_{r}}{\partial y},\varphi_{i}\right)\biggr{]}-\tilde{A}\,\bm{a}^{DNS}(t_{j})\biggr{\|}^{2}\,.$$

$$a^{DNS}_{i}(t_{j})=\biggl{(}\omega^{DNS}(t_{j}),\varphi_{i}\biggr{)}\,.$$

$$\dot{\mathbf{a}}=\mathbf{b}+(A+\tilde{A})\mathbf{a}+{\mathbf{a}}^{\top }B\mathbf{a},$$

$$\min_{\begin{subarray}{c}\tilde{A}\in\mathbbm{R}^{r\times r}\\[2.84544pt] \bm{a}^{\top}\tilde{A}\,\bm{a}\leq 0\end{subarray}}\,\ \sum_{j=1}^{M}\,\biggl{\|}\biggl{[}\left(\frac{\partial\omega_{R}}{\partial x}\frac{\partial\psi_{R}}{\partial y}-\frac{\partial\psi_{R}}{\partial x}\frac{\partial\omega_{R}}{\partial y},\varphi_{i}\right)-\left(\frac{\partial\omega_{r}}{\partial x}\frac{\partial\psi_{r}}{\partial y}-\frac{\partial\psi_{r}}{\partial x}\frac{\partial\omega_{r}}{\partial y},\varphi_{i}\right)\biggr{]}-\tilde{A}\,\bm{a}^{DNS}(t_{j})\biggr{\|}^{2}\,.$$

$$\bigl{\|}\overline{\psi^{DNS}({\bf x},\cdot)}-\overline{\psi^{ROM}({\bf x,\cdot})}\bigr{\|}^{2}_{L^{2}}\,\bigl{/}\,\bigl{\|}\overline{\psi^{DNS}({\bf x},\cdot)}\bigr{\|}^{2}_{L^{2}},$$

Full text

1. Introduction

Reduced order models (ROMs) for fluid dynamics have been abundantly investigated in recent decades as a way to reduce the computational cost of high resolution numerical schemes. The success of many ROM approaches has already been documented for various scientific and engineering applications, especially for flows that are governed by relatively few recurrent dominant spatial structures; see e.g. [3, 12, 16, 19, 21, 33, 38, 39, 45, 14] and references therein.

In this article, the recently proposed data-driven correction ROM (DDC-ROM) and its variants [53, 32, 22] are investigated in the numerical simulation of a quasi-geostrophic model of the double-gyre wind-driven ocean circulation. The DDC-ROMs fall into the category of hybrid projection/data-driven ROMs [10, 14, 34, 29, 20]. More specifically, in DDC-ROMs, the interactions among the resolved modes are the same as those in the standard Galerkin projection ROMs, while the interactions involving the unresolved modes are learned through a data-driven approach by fitting, e.g., a quadratic ansatz to the data that represents these missing interactions. In the following, we provide a brief derivation of the DDC-ROMs that builds on the standard projection ROMs, and refer to [53, 32] for more details.

To construct the standard projection ROM, we start with a general nonlinear system that has the following weak form111subject to possible further integration by parts for certain terms in $(\bm{f}(\bm{u})\,,\bm{v})$. in a suitable Hilbert space $\bm{X}$:

[TABLE]

where $\bm{f}$ is a general nonlinear function, $\bm{u}\in\bm{X}$ is the sought solution, and $(\cdot,\cdot)$ denotes the inner product on $\bm{X}$. Next, we use available data (snapshots) to construct orthonormal modes $\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{R}\}$, which represent the recurrent spatial structures, where $R$ is the rank of the snapshot matrix and typically $R=\mathcal{O}(10^{3})$ or even higher. Then, we choose the dominant modes $\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{r}\}$, typically with $r=\mathcal{O}(10)$, as ROM basis functions. The $r$-dimensional Galerkin ROM (G-ROM) of (1) is obtained by replacing $\bm{u}$ with a Galerkin truncation ${\bm{u}}_{r}=\sum_{j=1}^{r}a_{j}\,\bm{\varphi}_{j}$ and restricting $\bm{v}$ to the ROM subspace ${\bf X}^{r}:=\text{span}\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{r}\}$:

[TABLE]

In an offline stage, we construct the ROM, and in an online stage, we repeatedly use the G-ROM (2) with various parameters (if (1) is parameter dependent) and/or longer time intervals.

To construct the data-driven correction reduced order model (DDC-ROM) [53, 32, 22], we use an alternative approach: We start with a new Galerkin truncation, $\bm{u}_{R}=\sum_{j=1}^{R}a_{j}\,\bm{\varphi}_{j}$. We emphasize that, since $R=\mathcal{O}(10^{3})$ is the rank of the snapshot matrix, the new Galerkin truncation includes all the information in the available data (snapshots). Next, we replace $\bm{u}$ with $\bm{u}_{R}$ in (1) and project the resulting PDE onto ${\bf X}^{r}$:

[TABLE]

Since the ROM modes are orthonormal, \bigl{(}\overset{\bullet}{\bm{u}_{R}}\,,\bm{\varphi}_{i}\bigr{)}=\bigl{(}\overset{\bullet}{\bm{u}_{r}}\,,\bm{\varphi}_{i}\bigr{)},\,i=1,\ldots,r. Thus, (3) becomes

[TABLE]

which can be written as

[TABLE]

The last term on the right-hand side of (5) is a Correction term

[TABLE]

Thus, (5) can be written as G-ROM + Correction.

We emphasize that (5) is expected to be more accurate than the G-ROM, since the former is constructed from $\bm{u}_{R}$, whereas the latter is constructed from $\bm{u}_{r}$, where $r=\mathcal{O}(10)\ll R=\mathcal{O}(10^{3})$. Note that (5) is not yet a closed system in ${\bm{u}}_{r}$, since the Correction term involves $\bm{u}_{R}$, which lives in a higher-dimensional space than ${\bf X}^{r}$. Thus, to obtain from (5) an efficient $r$-dimensional ROM, we make the ansatz

[TABLE]

where ${\bf g}$ is a generic function (e.g., polynomial) whose coefficients/parameters still need to be determined. Once ${\bf g}$ is determined, the ROM (5) with the Correction term replaced by ${\bf g}$ yields the data-driven correction ROM (DDC-ROM):

[TABLE]

To determine the coefficients/parameters of the function ${\bf g}$ used in (8), we use data-driven modeling [5, 28, 37], i.e., we solve the following least squares problem:

[TABLE]

The numerical investigations in [53, 32, 22] show that the DDC-ROM (8) is significantly more accurate than the standard G-ROM in the numerical simulation of two test problems: (i) the 1D Burgers equation with a small diffusion coefficient $\nu=10^{-3}$; and (ii) a 2D flow past a circular cylinder at Reynolds numbers $Re=100,Re=500$, and $Re=1000$.

The main goal of this paper is to investigate the new DDC-ROM (8) in the numerical simulation of the quasi-geostrophic equations (QGE), which represent a significantly more difficult test case than the Burgers equation and the 2D flow past a circular cylinder considered in [53, 32]. Indeed, for the 2D flow past a circular cylinder with the Reynolds number $Re=1000$, 4the projection of the velocity field onto the first $8$ POD modes captures more than $99\%$ of the kinetic energy. In contrast, for the QGE in the parameter regime investigated below, a much broader range of spatial scales are actively involved in the time evolution of the turbulent fluid field. Indeed, the total amount of kinetic energy captured by the leading POD modes increases much slower for the QGE investigated here: it requires $16$ POD modes to capture $90\%$ of the kinetic energy, $37$ modes for $95\%$ of the kinetic energy, and $49$ modes for $96\%$ of the kinetic energy.

Furthermore, given the challenges posed by the QGE, we investigate two improvements to the DDC-ROM (8): First, we study the role of adding physical constraints to the DDC-ROM [32], in which the model for the Correction term in (7) satisfies the same type of physical constraints as those satisfied by the underlying equations; see (25). We also investigate whether modeling the commutation error, i.e., the error that appears as a result of interchanging spatial differentiation and ROM spatial filtering (e.g., projection) [22], improves the DDC-ROM accuracy. Finally, we study the DDC-ROM when it is trained on a time interval that is shorter than the time interval over which it is tested.

1.1. Connections to Previous Work

The DDC-ROM belongs to the class of ROM closure models, which model the effect of the truncated ROM modes (i.e., $\{\bm{\varphi}_{r+1},\ldots,\bm{\varphi}_{R}\}$) on the resolved ROM modes, (i.e., $\{\bm{\varphi}_{1},\ldots,\bm{\varphi}_{r}\}$). ROM closure models were first proposed in the pioneering work of Lumley and his collaborators [21] and are currently witnessing a dynamic development in several new directions, e.g., ROM spatial filtering, large eddy simulation (LES), and variational multiscale (VMS) [1, 2, 4, 52, 40], Mori-Zwanzig (MZ) formalism [36], nonlinear autoregression, moving averages with exogenous inputs (NARMAX) [9, 29], multilevel approaches and empirical model reduction (EMR) [25, 30, 24], data-adaptive harmonic decomposition and multilayer Stuart-Landau models (DAH-MSLM) [6, 23], and the parameterizing manifold (PM) approach rooted in the approximation theory of local invariant manifolds [8, 7], to name just a few. Probably the most dynamic development has been in using available data and machine learning techniques to develop ROM closure models [20, 31, 35, 43, 42, 54, 51].

The DDC-ROM is a hybrid projection/data-driven ROM, in which the standard Galerkin method is used to model the terms involving only the resolved modes and available data is used to model only the ROM closure term. This parsimonious/minimalistic data-driven approach is made possible by using ROM spatial filtering (i.e., ROM projection) and an LES/VMS framework to isolate the ROM closure term, which is then approximated by using data. The DDC-ROM minimalistic data-driven framework is similar in spirit to the NARMAX [9, 29] and PM [8, 7] ROM closure models, although they differ in the way the closure terms are handled. The DDC-ROM centers around ROM spatial filtering, whereas in the NARMAX approach the closure terms are modeled using nonlinear autoregression moving average with the resolved modes as exogenous inputs, and the PM approach parameterizes explicitly the unresolved modes in terms of the resolved modes.

The rest of the paper is organized as follows: In Section 2, we present briefly the QGE and the construction of the corresponding DDC-ROM. In Section 3, we assess the performance of the DDC-ROM using two metrics: the time-averaged streamfunction and the kinetic energy. Finally, in Section 4, we summarize our findings and outline future research directions.

2. Data-Driven Correction ROM (DDC-ROM)

In this section, we present the construction of the DDC-ROM for the QGE.

2.1. Quasi-geostrophic Equations (QGE)

In what follows, we use the quasi-geostrophic equations (QGE) as a mathematical model:

[TABLE]

where $\omega$ is the vorticity, $\psi$ is the streamfunction, $Re$ is the Reynolds number, and $Ro$ is the Rossby number. As a test problem for numerical investigation, we consider the QGE (10)–(11) with a symmetric double-gyre wind forcing given by

[TABLE]

The single-layer QGE (10)–(11) (also known as a barotropic vorticity equation (BVE)), are a popular mathematical model for forced-dissipative large scale ocean circulation. The idealized double-gyre wind forcing setting has been often used to understand the wind-driven circulation, e.g., the role of mesoscale eddies and their effect on the mean circulation. ROMs for the QGE (10)–(11) have been used, e.g., in [11, 13, 41, 47, 48].

The spatial domain of the QGE is $\Omega=[0,1]\times[0,2]$ and the time domain is $[0,80]$. We assume that $\psi$ and $\omega$ satisfy homogeneous Dirichlet boundary conditions:

[TABLE]

The QGE (10)–(11) can be cast in the general form of the nonlinear equation (1) by choosing

[TABLE]

We refrain from giving the precise formulation of the functional space here since this is tangential to the numerical study carried out below. The interested readers can consult for instance [17, Chapter 11], for the case of Navier-Stokes equations in the streamfunction-vorticity formulation. We will make precise how each of the terms in (14) is computed numerically once a set of POD basis functions for the vorticity is computed based on the direct numerical simulation (DNS) data obtained from a spectral code; see Sections 2.2 and 3.1.

2.2. Standard Galerkin ROM (G-ROM)

In our investigation, the ROM basis is obtained by using the proper orthogonal decomposition (POD)[21, 33]. We note, however, that other bases could be used for this purpose as well; examples include the dynamic mode decomposition [46], the principal interaction patterns [18, 26, 27], and the HIGAMod [38]. See also [49, 11, 50] for recent surveys and relationships/comparisons between different modal decomposition approaches.

We focus here mainly on the functional form of the G-ROM and defer details about the POD basis construction to Section 3. To this end, given an $r$-dimensional ROM subspace $X^{r}$ spanned by the first $r$ POD basis functions for the vorticity $\omega$,

[TABLE]

the $r$-dimensional G-ROM takes the form of (2) with $f$ therein given by (14). Recall that the streamfunction $\psi$ in (14) is related to the vorticity $\omega$ through the Poisson equation (11) subject to homogeneous Dirichlet boundary conditions.

To further reduce the G-ROM to an explicit ODE system, one option would be to replace $\psi$ in (14) by $-\Delta^{-1}\omega$, with $\Delta^{-1}$ being the inverse of the Laplacian of $\omega$ subject to the aforementioned boundary conditions. But since one important metric we adopt to assess the performance of the ROMs concerns the time average of $\psi$, we decide to keep $\psi$ explicit in the ROM formulation, although either way would lead to the same $r$-dimensional ODE system. For this reason, we also introduce a reduced set of $r$ basis functions for $\psi$, which are subordinate to the above POD basis functions for $\omega$ in (15) via $\phi_{i}(x,y)=-\Delta^{-1}\varphi_{i}(x,y)$, i.e., they solve the following Poisson equation:

[TABLE]

Note that while the POD basis $\{\varphi_{i}\}$ for the vorticity $\omega$ is an orthonormal basis under $L^{2}$-inner product, the basis $\{\phi_{i}\}$ for the streamfunction $\psi$ is not orthogonal. Given the G-ROM approximation $\omega_{r}=\sum_{i=1}^{r}a_{i}(t)\varphi_{i}(x,y)$ of $\omega$, the corresponding $\psi$ is approximated by $\psi_{r}=\sum_{i=1}^{r}a_{i}(t)\phi_{i}(x,y)$, which results from the ansatz $\psi_{r}=-\Delta^{-1}\omega_{r}$ and the above definition of the basis function $\phi_{i}$.

With the above notations, the $r$-dimensional G-ROM for the problem (10)–(13) is given by:

[TABLE]

where $(\cdot,\cdot)$ denotes the $L^{2}$ inner product over the spatial domain, and $i=1,\cdots,r$. Plugging the vorticity and streamfunction ROM approximations in (17), yields the Galerkin ROM (G-ROM):

[TABLE]

where \bm{a}\stackrel{{\scriptstyle def}}{{=}}\bigl{(}a_{j}(t)\bigr{)}_{j=1}^{r} is the vector of time-varying ROM coefficients. The G-ROM (18) can be written componentwise as follows: for $i=1,2,\cdots,r$

[TABLE]

where $b_{i}=Ro^{-1}(F,\varphi_{i}),\,A_{im}=Ro^{-1}\left(\frac{\partial\phi_{m}}{\partial x},\varphi_{i}\right)-Re^{-1}\left(\nabla\varphi_{m},\nabla\varphi_{i}\right),\,B_{imn}=-\left(J(\varphi_{m},\phi_{n}),\varphi_{i}\right).$

2.3. DDC-ROM

To construct the DDC-ROM (8) for the QGE, we adapt the general presentation in Section 1 to the QGE setting.

First, we note that the Correction term (7) takes the following form for the QGE:

[TABLE]

where $\omega_{R}({\bf x},t)=\sum_{i=1}^{R}a_{i}(t)\varphi_{i}({\bf x})$ and $\psi_{R}({\bf x},t)=\sum_{i=1}^{R}a_{i}(t)\phi_{i}({\bf x})$ are the $R$-dimensional ROM approximations of the vorticity and streamfunction in $X^{R}$, respectively, and ${\bf x}=(x,y)$. For clarity of presentation, in (20) we assume that the differentiation and the ROM projection commute (see, however, [22] for a detailed discussion of the commutation error). Thus, the linear terms in the QGE do not appear in the Correction (20). We emphasize that the Correction (20) is $R$-dimensional instead of $r$-dimensional, where $R\gg r$. Thus, to include the Correction (20) in the DDC-ROM, we first need to find an efficient, $r$-dimensional approximation of the Correction. To this end, we make the following linear ansatz: $\forall\,i=1,\ldots,r,$

[TABLE]

where the operator $\tilde{A}\in\mathbbm{R}^{r\times r}$ needs to be determined and \bigl{(}\tilde{A}\,\bm{a}\bigr{)}_{i} denotes the $i$-th component of the vector \bigl{(}\tilde{A}\,\bm{a}\bigr{)}. The ansatz (21) is chosen to resemble the right-hand side of the G-ROM (18); we note, however, that other ansatzes are possible [53, 32].

To compute the entries in the operator $\tilde{A}$ in (21), we use a data-driven approach. To this end, we adapt the least squares problem (9) to the QGE setting:

[TABLE]

In (22), $\bm{a}^{DNS}(t_{j})$ is the vector of ROM coefficients obtained from the DNS data, i.e., from the snapshots, at time instances $t_{j},\ j=1,\ldots,M$, which are obtained by projecting the corresponding snapshots $\omega^{DNS}(t_{j})=\sum_{k=1}^{R}a^{DNS}_{k}(t_{j})\,\varphi_{k}$ onto the POD basis functions $\varphi_{i}$ and using the orthogonality of the POD basis functions: $\forall\,i=1,\ldots,r,\ \forall\,j=1,\ldots,M,$

[TABLE]

The data-driven correction ROM (DDC-ROM) has the following form for the QGE:

[TABLE]

where the operators $\bm{b},A$, and $B$ are the G-ROM operators in (18) and the operator $\tilde{A}$ is the solution of the least squares problem (22).

3. Numerical Experiments

In this section, we investigate the DDC-ROM (24) in the numerical simulation of the QGE.

3.1. Computational Setting and Snapshot Generation

We investigate the QGE (10)–(11) with the symmetric double-gyre wind forcing given in (12) and homogeneous Dirichlet boundary conditions for $\psi$ and $\omega$ given in (13). The parameters are set to be $Re=450$ and $Ro=0.0036$.

For the DNS spatial discretization, we use a spectral method with a $257\times 513$ spatial resolution. Since both the vorticity and streamfunction have homogeneous boundary conditions, we approximate both functions with a tensor product Sine expansion in $x$ and $y$. For the DNS time discretization, we use an explicit Runge-Kutta method (Tanaka-Yamashita, an order $7$ method with an embedded order $6$ method for error control) and an error tolerance of $10^{-8}$ in time with adaptive time refinement and coarsening. We record the solution values every $10^{-2}$ simulation time units (starting at [math]) regardless of the current time step size so that the snapshots used in the POD are equally spaced. These spatial and temporal discretizations yield numerical results that are similar to the fine resolution numerical results obtained in [41, 44]. We follow [41, 44] and run the DNS between $[0,80]$. The time evolution of the spatially averaged kinetic energy in Fig. 1 in [41] shows that the flow converges to a statistically steady state, after a short transient interval that ends around $t=10$. We emphasize that even in the statistically steady state regime, the flow displays a high degree of variability. Thus, the numerical approximation of this statistically steady regime remains challenging for the low resolution ROMs that we investigate in this section.

3.2. ROM Construction

To generate the ROM basis (see Section 2), we follow [41, 44] and collect $701$ equally spaced snapshots of the vorticity, $\omega$, in the time interval $[T_{min},T_{max}]=[10,80]$ (on which the statistically steady state regime is attained) at equidistant time intervals. We also interpolate the DNS vorticity onto a uniform mesh with the resolution $257\times 513$ over the rectangle domain $\Omega=[0,1]\times[0,2]$, i.e., $h=\Delta x=\Delta y=1/256$. We then use the $701$ snapshots, form the correlation matrix $C$ for vorticity, and obtain the POD basis functions $\varphi_{i}$’s from the eigenvectors of $C$. Recall that the element $C_{ij}$ of $C$ is simply the $L^{2}$ inner product of the $i$-th and the $j$-th snapshots, i.e., $C_{ij}=\int_{\Omega}\omega_{i}\omega_{j}\mathrm{d}x\mathrm{d}y$ [41]. Throughout the article, the $L^{2}$ inner product over $\Omega$ is carried out by using the two-dimensional form of Simpson’s $1/3$ rule. Once the POD basis functions for the vorticity are generated, we solve the Poisson equation (16) to construct the streamfunction basis functions. For this purpose, we use a second order central difference (five-point stencil) spatial discretization of the Laplace operator.

In Fig. 1, we present the contour plots of selected streamfunction basis functions $\phi_{i}$’s and vorticity POD basis functions $\varphi_{i}$’s to give an idea about how the spatial scales are organized in these computed bases as the basis function index increases. Note that $\varphi_{i}$’s are much rougher than $\phi_{i}$’s, especially for the higher indices, whereas the roughness is smoothed out by the Laplacian when the $\phi_{i}$’s are computed according to (16). One source of roughness could be the uniform discretization mesh adopted here ($257\times 513$) when interpolating the DNS data, since there are steep gradients in the vorticity field as time evolves, both near the western boundary and within the domain. However, as we will illustrate in the next section, such roughness does not significantly degrade the DDC-ROM accuracy.

To construct the DDC-ROM, we need to solve the least squares problem (22), which can be ill-conditioned, especially when the training data is relatively short with respect to the number of coefficients that need to be learnt. To tackle this ill-conditioning issue, we use the truncated singular value decomposition (SVD) (see Step 6 of Algorithm 1 in [53]) with a tolerance that yields the most accurate results. Furthermore, to increase the computational efficiency of the DDC-ROM, we replace $(\omega_{R},\psi_{R})$ in (22) with $(\omega_{m},\psi_{m})$, where $1\leq m\leq R$ (for details, see Section 5.3 in [53]). Our numerical experiments suggest that $m=3r$ achieves a good balance between numerical accuracy and computational efficiency for the considered QGE model.

In our numerical investigation, in addition to the DDC-ROM, we also consider the physically-constrained DDC-ROM (CDDC-ROM) [32], which aims at improving the physical accuracy of the DDC-ROM. To construct the CDDC-ROM, we add physical constraints that require that the data-driven CDDC-ROM operators satisfy the same type of physical laws as those satisfied by the QGE. Specifically, we require that the CDDC-ROM’s Correction term’s linear component (i.e., the matrix $\tilde{A}$) should be dissipative. To implement these physical constraints, in the data-driven modeling step, we replace the unconstrained least squares problem (22) with a constrained least squares problem:

[TABLE]

Additionally, we also monitor the commutation error, which represents the effect of interchanging ROM spatial filtering and differentiation [22]. For this test problem, modeling the commutation error does not significantly change the DDC-ROM and CDDC-ROM results, suggesting that the commutation error does not play a significant role in the ROM construction. Thus, for clarity of presentation, we do not include the commutation error in the DDC-ROM and CDDC-ROM results.

In the online stage, for all the ROMs we utilize the fourth order Runge-Kutta scheme (RK4) for the temporal discretization. To ensure the numerical stability of the time discretization, we choose a time step size $\Delta t=0.001$. We store ROM data every ten time steps to match the DNS sampling rate. We use the DNS snapshot at $t=10$ to initialize the ROMs.

3.3. Numerical Results

In this section, we present numerical results for the G-ROM, DDC-ROM, and CDDC-ROM. As benchmark for our numerical investigation, we use the DNS results.

3.3.1. Kinetic energy

In this section, we assess the performance of the ROMs using the DNS kinetic energy as a metric. As pointed out in Section 1, due to the involvement of a broad range of active spatial scales, it would be too demanding to require any ROMs to reproduce the statistics of the DNS kinetic energy or any other reasonable observables when the dimension is too low, at least within the POD basis framework adopted here. Thus, we confine ourselves instead to a much less ambitious goal of reproducing the range of oscillations presented in the DNS kinetic energy. The assessment at a more quantitative level will be carried out in the next section for another metric.

Recall that the velocity field $(u(\bm{x},\cdot),v(\bm{x},\cdot))$ used in the computation of the kinetic energy E(t)=\frac{1}{2}\int_{\Omega}\big{(}u^{2}(\bm{x},t)+v^{2}(\bm{x},t)\big{)}d\bm{x} is related to the streamfunction according to $(u,v)=(\partial_{y}\psi,-\partial_{x}\psi)$. The first-order spatial derivatives are calculated using a $4$-th order accurate central difference scheme. The kinetic energy itself is computed using the two-dimensional form of Simpson’s $1/3$ rule.

In Fig. 2, for three different $r$ values ($r=10,15$, and $40$), we plot the time evolution of the ROM kinetic energy. For $r=10$, the G-ROM kinetic energy takes off very quickly and stabilizes at a level around $8\times 10^{4}$, which is roughly 200 times higher than the DNS kinetic energy on average. In contrast, both the DDC-ROM and CDDC-ROM successfully stabilize the G-ROM, and produce kinetic energies almost of the same order of magnitude as the DNS kinetic energy (although there is some overdamping in the CDDC-ROM result due to the physical constraint $\bm{a}^{\top}\tilde{A}\,\bm{a}\leq 0$; see (25)).

For $r=15$, the G-ROM kinetic energy is within good range at the beginning of the simulation, but increases to an unphysical value around $t=40$ and eventually stabilizes and oscillates around $2\times 10^{4}$. In contrast, the DDC-ROM and CDDC-ROM kinetic energies are both within the good range, and between the two, the CDDC-ROM performs better in reproducing the peaks and the corresponding frequency of the peaks.

Finally, for $r=40$, the G-ROM and DDC-ROM perform similarly. The DDC-ROM kinetic energy is closer to the DNS kinetic energy over certain time windows (e.g., $[50,60]$), whereas the CDDC-ROM kinetic energy is somewhat lower than the DNS kinetic energy.

The above results suggest that both the DDC-ROM and CDDC-ROM can successfully stabilize a severely truncated G-ROM. For the chosen criterion, the advantage of the DDC-ROM over the G-ROM is clearly visible for all ROM dimensions $r$ between $10$ and $30$. The CDDC-ROM can produce results comparable or even better than DDC-ROM for $r$ between $10$ and $20$. For even higher dimensions, the CDDC-ROM tends to overdamp the kinetic energy. This is plausible, since the physical constraint $\bm{a}^{\top}\tilde{A}\,\bm{a}\leq 0$ in the estimation of the matrix $\tilde{A}$ for the CDDC-ROM aims to enhance the stability of the ROM, but does not also guarantee improved accuracy compared to the DDC-ROM.

3.3.2. Relative errors for the time-averaged streamfunction

In this section, we assess the ROM performance at a more quantitative level using the ROM time-averaged streamfunction over the aforementioned time interval, $[10,80]$. It is known that the time-averaged streamfunction displays a four-gyre structure [15] even though a double-gyre wind forcing is employed; cf. (12). The metric that we use is the following relative error:

[TABLE]

where $\overline{(\cdot)}$ represents time average over $[10,80]$, and ${\bf x}=(x,y)$.

In Table 1, we list this relative error for each of the ROMs as the ROM dimension $r$ increases. For small $r$ values (i.e., $5\leq r\leq 20$), the CDDC-ROM is the most accurate. Indeed, for $r=5$, the CDDC-ROM is the only ROM that yields a stable approximation: all other ROMs with $r=5$ experience exponential blowup with the given timestep. Furthermore, for $r=10$ and $r=15$, the CDDC-ROM error is at most half of the DDC-ROM error. Finally, for $5\leq r\leq 50$, the G-ROM error is one to three orders of magnitude larger than the CDDC-ROM error. These results are also supported by the plots in Fig. 3, which display the time-average of the streamfunction $\psi$ over the time interval $[10,80]$ for DNS, G-ROM ($r=10$), DDC-ROM ($r=10$), and CDDC-ROM ($r=10$). These plots clearly show that the DDC-ROM and CDDC-ROM are able to capture the correct four-gyre structure, whereas the G-ROM fails drastically at this low ROM dimension.

For large $r$ values (i.e., $25\leq r\leq 50$), the DDC-ROM results in Table 1 are the most accurate. Indeed, at $r=25$ the CDDC-ROM error starts to increase, whereas the DDC-ROM error generally decreases. The G-ROM error also continues to decrease, but it is always larger than the DDC-ROM error.

We conclude that the CDDC-ROM is the most accurate for small $r$ values and the DDC-ROM is the most accurate for large $r$ values. These results suggest that adding physical constraints to the DDC-ROM is beneficial in the highly truncated cases (i.e., for small $r$ values), but the benefit brought by the physical constraints diminishes as $r$ is further increased, and can even produce less accurate results than DDC-ROM due to overdamping, as pointed out in Section 3.3.1. We emphasize, however, that we are using the linear ansatz to construct the DDC-ROM; further numerical investigations are needed to determine the role of physical constraints when the DDC-ROM is built with higher-order (e.g., quadratic) ansatz [32]. We also note that the G-ROM is consistently less accurate for all the $r$ values. Finally, we note that, as $r$ increases, the errors for all the ROMs reach a plateau instead of converging to zero. We believe that this behavior is due to the roughness present in the vorticity basis functions, especially for the higher indices (see Fig. 1).

3.3.3. Shorter training time interval

In this subsection, we consider the situation when the ROMs are trained on a time interval that is shorter than the time interval over which the ROMs are tested. Specifically, we only use snapshots in the time interval $[10,t_{p}^{\ast}]$, sampled every 0.1 time units as before, to generate the ROM basis and construct the ROM operators $A,B,$ and $\tilde{A}$. This leads to a total number of $(t_{p}^{\ast}-10)/0.1+1$ snapshots. We investigate two different cases: (I) $t_{p}^{\ast}=45$ and (II) $t_{p}^{\ast}=35$.

In Table 2, we list the relative errors associated with the ROMs for the time-averaged streamfunction defined in (26). The results show that the DDC-ROM is significantly more accurate than the G-ROM, especially for small $r$ values. This is in line with the results in Table 1. Furrthermore, we note that, as expected, the shorter the time interval $[10,t_{p}^{\ast}]$ is, the larger the ROM relative error.

4. Conclusions

We enhanced the standard Galerkin ROM (G-ROM) for the quasi-geostrophic equations (QGE) with an additional term derived from available data and a least squares optimization procedure. These ideas are based on previous work and are usually referred to as data-driven correction ROMs (DDC-ROMs) and constrained data-driven correction ROMs (CDDC-ROMs). The latter incorporate a negative semidefiniteness constraint into the optimization problem to preserve a fundamental property of the linear operator in the DDC-ROMs.

The QGE are challenging equations that exhibit complex spatiotemporal behavior: we were able to significantly improve the G-ROM performance by adding an additional term (derived by optimization) to the linear component of the G-ROM. For a ROM with $10$ POD modes, the DDC-ROM lowered the error in the mean streamfunction (compared to the G-ROM) by a factor of about $600$; similarly, the CDDC-ROM lowered the error by a factor of about $2000$.

In the future, we plan on investigating whether using a higher-order (e.g., quadratic) ansatz in the construction of the DDC-ROM and CDDC-ROM yields more accurate results than using a linear ansatz (i.e., the approach utilized in this paper). We also plan to study parameter sensitivity (on both the Reynolds and Rossby numbers) and examine the possibility of constructing a sequence of ROMs that work across a wide range of either value.

Acknowledgments

This work is partially supported by NSF DMS-1821145 (CM,TI), NSF DMS-1616450 (HL), and NSF DMS-1450327 (DRW).