Continuous Data Assimilation Reduced Order Models of Fluid Flow

TL;DR

This paper introduces a novel continuous data assimilation reduced order model (DA-ROM) for fluid flow simulation, improving long-term accuracy by integrating measurement data and adaptive nudging strategies, with proven exponential convergence.

Contribution

The paper combines continuous data assimilation with reduced order models, providing a new approach that enhances long-term accuracy and stability in fluid flow simulations.

Findings

DA-ROM converges exponentially fast to the true solution with proper nudging.

Adaptive nudging improves long-term accuracy of ROMs.

Numerical tests confirm the effectiveness of DA-ROM with adaptive strategies.

Abstract

We propose, analyze, and test a novel continuous data assimilation reduced order model (DA-ROM) for simulating incompressible flows. While ROMs have a long history of success on certain problems with recurring dominant structures, they tend to lose accuracy on more complicated problems and over longer time intervals. Meanwhile, continuous data assimilation (DA) has recently been used to improve accuracy and, in particular, long time accuracy in fluid simulations by incorporating measurement data into the simulation. This paper synthesizes these two ideas, in an attempt to address inaccuracies in ROM by applying DA, especially over long time intervals and when only inaccurate snapshots are available. We prove that with a properly chosen nudging parameter, the proposed DA-ROM algorithm converges exponentially fast in time to the true solution, up to discretization and ROM truncation errors. Finally, we propose a strategy for nudging adaptively in time, by adjusting dissipation arising from the nudging term to better match true solution energy. Numerical tests confirm all results, and show that the DA-ROM strategy with adaptive nudging can be highly effective at providing long time accuracy in ROMs.

Tables

Table 1. Table 1: DA-ROM rates of convergence with respect to the ROM truncation.

No. modes	${\left({\sum }_{j=r+1}^d{\lambda }_j\left(1+{\Vert \nabla {\phi }_j\Vert }^2\right)\right)}^{1/2}$	Error	Rate
8	2.218e+2	4.980e-2	–
10	1.077e+2	4.850e-2	1.74
12	9.246e+1	3.046e-2	2.51
14	7.680e+1	1.793e-2	1.70
16	4.590e+1	1.360e-2	1.36
18	3.334e+1	9.498e-3	1.12
20	2.601e+1	6.974e-3	1.24

Equations

$$\mu I_H(u_r-u_{obs}),$$

$$u_t+u\cdot \nabla u+\nabla p-\nu \Delta u=f,\text{and }\nabla \cdot u=0,\text{in}\Omega \times (0,T]$$

$$u=0,\text{on}\partial \Omega \times (0,T],\text{and }u(x,0)=u_0(x),\text{in}\Omega .$$

$$V_h:=\{v_h\in X_h:(\nabla \cdot v_h,q_h)=0\forall q_h\in Q_h\}\subset X.$$

$$\Vert \varphi \Vert \le C_P\Vert \nabla \varphi \Vert \forall \varphi \in X.$$

$$b(w,u,v)=(w\cdot \nabla u,v)\forall u,v,w\in X,$$

$$b^{\ast }(w,u,v):=\frac{1}{2}(w\cdot \nabla u,v)-\frac{1}{2}(w\cdot \nabla v,u)\forall u,v,w\in X.$$

$$|b^{\ast }(u,v,w)|$$

$$r{a}_{n+1}\le a_n+B,$$

$$a_{n+1}\le a_0{\left(\frac{1}{r}\right)}^{n+1}+\frac{B}{r-1}.$$

$$u_r(x,t)=\sum _{j=1}^r{a}_j(t){\phi }_j(x),$$

$$(u_{r,t},{\phi }_i)+\nu (\nabla u_r,\nabla {\phi }_i)+b^{\ast }(u_r,u_r,{\phi }_i)=(f,{\phi }_i).$$

$$(P_r(v),v_r)=(v,v_r)\forall v_r\in X_r.$$

$$\Vert \nabla \phi \Vert \le |||{\mathbb{S}}_R||{|}_2^{1/2}\Vert \phi \Vert \forall \phi \in X_r,$$

$$\Vert u^M-u_h^M{\Vert }^2+\nu h^2\Delta t\sum _{n=1}^M\Vert \nabla (u^n-u_h^n){\Vert }^2$$

$$\sum _{n=1}^M\Vert u^n-P_r(u^n){\Vert }^2$$

$$\sum _{n=1}^M\Vert \nabla (u^n-P_r(u^n)){\Vert }^2$$

$$\displaystyle\hskip 113.81102pt+\sum_{j=r+1}^{d}\|\nabla{\varphi}_{j}\|^{2}\lambda_{j}\bigg{)}.$$

$$\underset{n}{\max }\Vert u^n-P_r(u^n){\Vert }^2$$

$$\underset{n}{\max }\Vert \nabla (u^n-P_r(u^n)){\Vert }^2$$

$$\displaystyle\hskip 113.81102pt+\sum_{j=r+1}^{d}\|\nabla{\varphi}_{j}\|^{2}\lambda_{j}\bigg{)}.$$

$$\Vert u^M-u_h^M{\Vert }^2+h^2\Vert \nabla (u^n-u_h^n){\Vert }^2\le C(\nu ,p)(h^{2k+2}+\Delta t^4),$$

$$\Vert I_H(w)-w\Vert$$

$$\Vert I_H(w)\Vert$$

$$\frac{1}{\Delta t}(u_r^{n+1}-u_r^n,v_r)+b^{\ast }(u_r^{n+1},$$

$$+\mu (I_H(u_r^{n+1}-u(t^{n+1})),I_H{v}_r)=(f^{n+1},v_r),$$

$$\Vert u_r^M{\Vert }^2\le \Vert u_r^0{\Vert }^2{\left(\frac{1}{1+\lambda \Delta t}\right)}^M+C{\lambda }^{-1}({\nu }^{-1}{F}^2+\mu U^2):=C_{data},$$

$$u\in L^{\infty }(0,T;H^1(\Omega ))\cap H^1(0,T;H^{k+1}(\Omega ))\cap H^2(0,T;H^1(\Omega )),$$

$$p\in L^2(0,T;H^{k+1}(\Omega )),$$

$$f\in L^2(0,T;L^2(\Omega )).$$

$${\alpha }_1$$

$${\alpha }_2$$

$$\Vert u^{n+1}$$

$$\displaystyle+C\lambda^{-1}\bigg{\{}{\Delta t}^{2}+\nu^{-1}h^{2k}+\beta_{1}C^{2}_{I}\mu\bigg{(}h^{2k+2}+{\Delta t}^{4}+\sum_{j=r+1}^{d}\lambda_{j}\bigg{)}$$

$$\displaystyle+(\nu^{-1}M^{2}+\nu^{-1}M^{2}|||\mathbb{S}_{R}|||_{2})\bigg{(}h^{2k}+{\Delta t}^{4}+\sum_{j=r+1}^{d}\|\nabla\varphi_{j}\|^{2}\lambda_{j}\bigg{)}\bigg{\}},$$

Full text

1 Introduction

Reduced order models (ROMs) for fluids dominated by relatively few recurrent spatial structures are generally built as follows [14, 15, 33]: (i) postulate a collection of snapshots, either from numerical experiments or from physical data; (ii) from those snapshots, select a small number (e.g., $10$) of ROM basis functions; (iii) project the equations of motion into this basis; and (iv) advance the velocity in time to interrogate flows different from the one generating the snapshots. ROMs have been explored for decades [15]. When successful, ROMs can decrease the computational cost of a brute force, direct numerical simulation (DNS) by orders of magnitude.

One of the main roadblocks for ROMs of realistic flows is their lack of accuracy, e.g., in complex problems, for long time intervals, or when a low-dimensional ROM basis is used. To increase the ROM accuracy in practical applications, several approaches are currently used. We list some of these below:

(i) Closure Modeling:  To model the effect of the discarded ROM modes, a Correction term is generally added to the standard ROM [3, 29, 32, 37]. Given the drastic truncation used in ROMs for realistic flows, the Correction term is essential for accuracy.

(ii) Numerical Stabilization:  To eliminate/alleviate the spurious numerical oscillations generated when ROMs are used for convection-dominated flows, numerical stabilization techniques can be used [7, 13, 37].

(iii) Data-Driven Modeling:  Recently, available numerical or experimental data have been used to construct ROM operators [30] or to determine the unknown coefficients in classical ROM operators [3, 11, 32].

(iv) Improved Basis:  Another approach for increasing the ROM accuracy in practical applications is the construction of an improved (more accurate) ROM basis that better captures the behavior of the underlying system [2, 35, 36, 41].

(v) Physical Accuracy:  To develop physically sound ROMs, recent effort has been directed at ensuring that the ROMs satisfy the same physical balances/conservation laws as those satisfied by the equations of motion [28, 27].

In this paper, we propose a new approach to increase the ROM accuracy. Specifically, we use data assimilation (DA) to develop a novel DA-ROM. In weather modeling, climate science, and hydrological and environmental forecasting, DA has been used for decades to incorporate observational data in simulations, in order to increase the accuracy of solutions and to obtain better estimates of initial conditions [20]. In this paper, we use DA to improve the ROM accuracy. Specifically, we add to the standard ROM a feedback control term of the form

[TABLE]

which nudges the ROM approximation ($u_{r}$) towards the reference solution ($u_{obs}$) corresponding to the observed data. In (1.1), $I_{H}$ is an interpolation operator onto a coarser mesh of size $H$ and $\mu>0$ is a nudging parameter. Equation (1.1) allows the simple implementation of DA into existing ROM codes. The nudging term (1.1) increases the accuracy of the new DA-ROM by utilizing the available low-resolution data, without the need to increase the number of ROM basis functions. This type of data assimilation has recently become popular due to a seminal paper of Azouani et al. [1], and since then it has been used to improve solutions to many different types of evolutionary systems [4, 5, 8, 9, 10, 18, 22, 26].

We emphasize that our new DA-ROM is different from other uses of DA for ROMs, e.g., [6, 19, 25, 39]: In the latter the authors use ROMs to speed up classical DA algorithms (e.g., 4D-VAR), whereas in the DA-ROM proposed in this paper, we use DA to improve the ROM accuracy.

The rest of the paper is organized as follows: In Section 2, we introduce some notation and preliminaries necessary for our analysis. In Section 3, we construct the new DA-ROM and perform a careful error analysis. In Section 4, we perform a numerical investigation of the new DA-ROM in the numerical simulation of a 2D flow past a circular cylinder and discuss implementation of the DA-ROM algorithm with an adaptive nudging parameter, which can be used to further improve the accuracy of solutions. Finally, in Section 5, we draw conclusions and outline future research directions.

2 Notation and Preliminaries

Let $\Omega\subset\mathbb{R}^{d}$, $d$=2 or 3, be a bounded open domain. The $L^{2}(\Omega)$ norm and inner product will be denoted by $\|\cdot\|$ and $(\cdot,\cdot)$, respectively, and all other norms will be appropriately labeled with subscripts.

We consider the Navier-Stokes equations (NSE) with no-slip boundary conditions:

[TABLE]

Here $u$ is the velocity, $f=f(x,t)$ is the known body force, $p$ is the pressure, and $\nu$ is the kinematic viscosity.

We denote the natural velocity space by $X=H^{1}_{0}(\Omega)$ and pressure space by $Q=L_{0}^{2}(\Omega)$ and by $(X_{h},Q_{h})\subset(X,Q)$, corresponding inf-sup stable finite element spaces. Additionally, we define the discretely divergence-free space $V_{h}$ as

[TABLE]

The Poincaré inequality will be used throughout this paper: there exists a constant $C_{P}$ depending only on $\Omega$ such that

[TABLE]

We define the trilinear form

[TABLE]

and the explicitly skew-symmetric trilinear form given by

[TABLE]

An important property of the $b^{\ast}$ operator is that $b^{\ast}(u,v,v)=0$ for $u,v\in X$. We will utilize the following bounds on the operator $b^{*}$ [24].

Lemma 2.2**.**

There exists a constant $M>0$ dependent only on $\Omega$ satisfying

[TABLE]

for all $u,v,w\in X$ for which the norms on the right hand sides are finite.

The following lemma is proven in [23], and is useful in our analysis.

Lemma 2.3**.**

Suppose constants $r$ and $B$ satisfy $r>1$, $B\geq 0$. Then if the sequence of real numbers $\{a_{n}\}$ satisfies

[TABLE]

we have that

[TABLE]

2.1 ROM preliminaries

Let $\{u_{h}^{1},...,u_{h}^{M}\}$ be snapshots of FE solutions at $M$ different time instances. The proper orthogonal decomposition seeks a low-dimensional basis that approximates these snapshots optimally with respect to a certain norm; in this paper, we use the $L^{2}$ norm. This minimization can be set up as an eigenvalue problem $YY^{T}M_{h}\varphi_{j}=\lambda_{j}\varphi_{j}$, $j=1,...,N_{h}$. where $N_{h}$ is the dimension of the finite element space. The eigenvalues are real and non-negative, so they can be ordered as $\lambda_{1}\geq...\geq\lambda_{d}\geq\lambda_{d+1}=...=\lambda_{N_{h}}=0$, where $d$ is the rank of the snapshot matrix. We take the ROM space to be $X_{r}:=\text{span}\{\varphi_{i}\}_{i=1}^{r}$, and note that $X_{r}\subset V_{h}$. The ROM approximation of the velocity is defined as

[TABLE]

where the coefficients $a_{j}(t)$ are determined by solving the Galerkin ROM:

[TABLE]

We define the $L^{2}$ ROM projection $P_{r}:L^{2}\to X_{r}$ by: for all $v\in L^{2}(\Omega)$, $P_{r}(v)$ is the unique element of $X_{r}$ such that

[TABLE]

In addition, the following inverse inequality holds for our ROM basis [21].

Lemma 2.6** (POD inverse estimate).**

[TABLE]

where $|||\mathbb{S}_{R}|||$ is the matrix 2-norm of the ROM stiffness matrix, as in Lemma 3.1 of [16].

In order to establish an error estimate for the ROM projection, we first make the following assumption on the finite element error:

Assumption 2.8**.**

Let $C(\nu,p)$ denote a constant which is dependent upon the viscosity and pressure. We assume that the finite element error $u_{h}$ satisfies the following error estimate

[TABLE]

Remark 2.10**.**

Error estimates of this form have been proven for varying amounts of regularity on the continuous solution $u$ and $p$. Some examples include the scheme used in the numerical experiments in Section 4.

Using Assumption 2.8 the following error estimates for the ROM projection can be proven [16]:

Lemma 2.11**.**

The $L^{2}$ ROM projection of $u^{n}$ satisfies the following error estimates:

[TABLE]

We then make the following assumption similar to that made in [16]:

Assumption 2.14**.**

The $L^{2}$ ROM projection of $u^{n}$ satisfies the following error estimates:

[TABLE]

Remark 2.17**.**

If we assumed in Assumption 2.8 that the finite element error satisfies

[TABLE]

then the bound in Assumption 2.14 would hold. Error estimates of this form have been proven for varying amounts of regularity on the continuous solution $u$ and $p$. Some examples include the incremental pressure correction schemes in [12] and chapter 7 of [31].

2.2 Data assimilation preliminaries

We consider $I_{H}$ to be an interpolation operator that satisfies: For a given mesh $\tau_{H}(\Omega)$ with $H\leq 1$,

[TABLE]

for any $w\in H^{1}(\Omega)$. For example, this holds for the $L^{2}$ projection onto piecewise constants, and the Scott-Zhang interpolant. For the (unknown) true solution $u$, $I_{H}(u)$ represents an approximation of what is observed of the true solution. We assume in this paper that $I_{H}(u)$ can be observed at any time. We remark that in the nudging used in this paper, we use the interpolant of the test function as well, as suggested in [34]. This does not affect the convergence analysis; some extra terms arise that can be handled without difficulty, but it allows for unconditional stability.

3 Error Analysis

For simplicity of exposition, our analysis considers a first order DA-ROM algorithm, which takes the following form: Find $u_{r}^{n+1}\in X_{r}$ such that for all $v_{r}\in X_{r}$,

[TABLE]

for $n=1,2,...,M$, with $v_{0}=P_{r}(u_{0})$, and where $\mu\geq 0$ is the nudging parameter. Extension to other time stepping methods is possible, and, for example, extension to BDF2 can be done following the usual techniques [23]. All of our numerical tests use the analogous BDF2 algorithm.

We first prove a stability estimate for the DA ROM algorithm.

Lemma 3.2**.**

The solutions to (3.1) satisfy for all $M>1$,

[TABLE]

where $F:=\|f\|_{L^{\infty}(0,\infty;H^{-1})}$, $U:=\|u\|_{L^{\infty}(0,\infty;L^{2})}$, and $\lambda=\nu C_{P}^{-2}$.

Proof.

This result follows as in [23] by letting $v_{r}=u_{r}^{n+1}$ in (3.1) and using Cauchy-Schwarz and Young’s inequalities. Additionally, the non-negative DA term $\|I_{H}u_{r}^{n+1}\|^{2}$ can be dropped from the left after bounding the right hand side. ∎

To analyze rates of convergence of the approximation we make the following regularity assumptions on the NSE [24]:

Assumption 3.3**.**

We assume that the solution of the NSE satisfies

[TABLE]

We next prove that solutions to (3.1) converge to the true solution exponentially fast, up to discretization and ROM projection error.

Theorem 3.4**.**

Define

[TABLE]

which have parameters $\mu$, $H$, $\beta_{i}>0$, $i=1,2$ that are chosen so that $\alpha_{i}>0$, $i=1,2$. Then under the regularity assumptions of Assumption 3.3, we have that

[TABLE]

where $\lambda=\min\{\alpha_{1}C_{P}^{-2},\alpha_{2}\}$.

Remark 3.6**.**

The ${\Delta t}^{2}$ term that shows up on the right hand side of (3.5) is a result of the first order time stepping in Algorithm 3.1. If we instead used a second order approximation, like BDF2, then this term would be replaced by ${\Delta t}^{4}$.

Remark 3.7**.**

If $\beta_{1},\beta_{2}$ are chosen to be $1/2$, the condition $\alpha_{1}>0$ reduces to $\nu-C\mu H^{2}>0$, which is the same condition found in [23] and references therein, for a relationship between the nudging parameter, viscosity, and coarse mesh width. Choosing $\beta_{1},\beta_{2}$ larger can allow one to choose the coarse mesh width $H$ larger (and thus require less observational data) while still satisfying $\alpha_{i}>0$, $i=1,2$. However, there is a trade-off because $\beta_{1}$ appears on the right hand side of equation (3.5): As $\beta_{1}$ increases, so does the bound on the DA-ROM error.

Proof.

The NSE (true) solution satisfies

[TABLE]

Note that we can write the time derivative term above as $C{\Delta t}u_{tt}(t^{*})$ for some $t^{*}\in(t^{n},t^{n+1})$ [23]. Subtracting (3.1) from (3.8) and letting $e^{n}:=u_{r}^{n}-u^{n}$, we obtain

[TABLE]

Decompose the error as a part inside the ROM space and one outside by adding and subtracting the $L^{2}$ projection of $u^{n}$ into the ROM space, $P_{r}(u^{n})$ (see (2.5)):

[TABLE]

Letting $v_{r}=\phi_{r}^{n+1}$ in (3.9), we note that since $\phi_{r}^{n+1}\in X_{r}\subset V_{h}$, for any $q_{h}\in Q_{h}$,

[TABLE]

Adding and subtracting $\phi_{r}^{n+1}$ to both components of the nudging term we have

[TABLE]

Using the polarization identity, the fact that $(\eta^{n+1}-\eta^{n},\phi^{n+1}_{r})=0$ (by the definition of the $L^{2}$ projection), and dropping the nonnegative term $\frac{1}{2{\Delta t}}\|\phi_{r}^{n+1}-\phi_{r}^{n}\|^{2}$ on the left hand side, we have

[TABLE]

By Poincaré, Cauchy Schwarz, and Young’s inequalities, we bound the first two terms on the right hand side and the pressure term,

[TABLE]

The first two nonlinear terms are now bounded similarly to those in [28] using Cauchy-Schwarz and Young’s inequalities, and the first inequality from Lemma 2.2:

[TABLE]

How we treat the third nonlinear term is the key difference in the proof from standard schemes (see chapter 9 of [24]). Due to the added dissipation from the DA term on the left-hand side of (3.12), we will be able to hide the term containing $\phi^{n+1}_{r}$, rather than invoking a discrete Gronwall’s inequality. Thus, for this term we use the second inequality from Lemma 2.2 and the ROM inverse inequality (2.7) to obtain

[TABLE]

The first nudging terms on the right hand side of (3.12) are bounded using (2.19), Cauchy Schwarz, and Young’s inequality

[TABLE]

The second nudging term is bounded using Cauchy Schwarz and Young’s inequality, and (2.18), yielding

[TABLE]

Now letting $c_{i}=\frac{1}{12}$, $i=1,2,..,6$, combining terms, and recalling our definition of $\alpha_{1}$ and $\alpha_{2}$ given in the statement of the theorem, (3.12) becomes

[TABLE]

where $C$ is a generic constant which is independent of $\nu,p,u,T,H,C_{I}$. Next, we bound the fourth term on the right hand side further using the ROM inverse inequality (2.7), and the stability result from Lemma 3.2

[TABLE]

Now applying Lemma 2.11, using our regularity assumptions, and taking $\lambda:=\min\{C_{P}^{-2}\alpha_{1},\alpha_{2}\}$ in (3.19), it then follows that

[TABLE]

Finally, by Lemma 2.3, we obtain

[TABLE]

The triangle inequality completes the proof.

∎

4 Numerical Experiments

In this section, we perform a numerical investigation of the new DA-ROM. In Section 4.1, we illustrate the theoretical scalings proved in Section 3. In Section 4.2, we investigate the numerical accuracy of the new DA-ROM. In Section 4.3, we investigate the new DA-ROM when inaccurate snapshots are used in its construction. Finally, in Section 4.4, we propose and investigate an adaptive nudging procedure.

We consider Algorithm 3.1 (except here with BDF2) applied to 2D channel flow past a cylinder [38], with Reynolds number $Re$=500. The domain is the rectangular channel [0, 2.2]$\times$[0, 0.41], with a cylinder centered at $(0.2,0.2)$ and radius $0.05,$ see Figure 1. There is no external forcing ($f=0$), no-slip boundary conditions are prescribed for the walls and the cylinder, and an inflow profile is given by

[TABLE]

We take $\nu=0.0002$, and enforce the zero-traction boundary condition with the usual ‘do-nothing’ condition at the outflow.

The DNS is run to t=15 with the usual BDF2-FEM discretization [24] using $(P_{2},P_{1}^{disc})$ Scott-Vogelius elements on a barycenter refined Delaunay mesh that provided 103K velocity dof, a time step of $\Delta t=0.002$, and with the simulations starting from rest ($v_{h}^{0}=v_{h}^{1}=0$). Lift and drag calculations were performed for the computed solution and compared to the literature [38, 40], which verified the accuracy of the DNS. We used the snapshots from t=5 to t=6 to generate the ROM modes.

For the lift and drag calculations, we used the definitions

[TABLE]

where $u_{t_{S}}$ is tangential velocity, $S$ is the cylinder, and $n=\langle n_{x},n_{y}\rangle$ the outward unit normal to the domain. For the calculations, we used the global integral formulation from [17].

The coarse mesh for DA is constructed using the intersection of a uniform rectangular mesh with the domain. We take $H$ to be the width of each rectangle, and use $H=\frac{2.2}{20}$ (400 measurement locations) in our tests. Figure 2 shows in red a 35K dof mesh and associated $H=\frac{2.2}{8}$ coarse mesh in black.

For the DA-ROM computations, we start from zero initial conditions $v_{h}^{1}=v_{h}^{0}=0$, use the same spatial and temporal discretization parameters as the DNS, and start assimilation with the $t=5$ DNS solution (i.e., time 0 for DA-ROM corresponds to $t=5$ for the DNS).

4.1 Rates of Convergence

In this section, we illustrate numerically the rates of convergence in Section 3. Theorem 3.4 gave a DA-ROM error estimate that depends on the ROM eigenvalues and eigenvectors, for sufficiently large $n$ and assumptions on $\mu$ and $H$:

[TABLE]

where $(\lambda_{j},\varphi_{j})$ are the eigenpairs of the ROM eigenvalue problem described in Section 2.1. Table 1 illustrates the dependence of the error bound on the dimension of the DA-ROM space, $r$. Taking $T$=1, $\mu=100$, $H=\frac{2.2}{20}$, $Re$=500, we run the ROM with varying $r$ and calculate the $L^{2}$ spatial error at the last time step. We also calculate the quantity in the error estimate corresponding to the eigenvalues and eigenvectors (i.e., $\sum_{j=r+1}^{d}\lambda_{j}(1+\|\nabla\varphi_{j}\|^{2}))^{1/2}$), and use this and the error to calculate the corresponding convergence rate with respect to increasing $r$. From the theorem, we expect a rate of 1, and our results are consistent with this rate.

4.2 Numerical Accuracy

In this section, we investigate the numerical accuracy of the new DA-ROM. Specifically, we compare the performance of the DA-ROM to that of the standard ROM ($\mu=0$) and the DNS solution in predicting energy and drag (lift is accurate in all of our tests, so we omit it here). We run to t=10, and run tests with both $N=8$ and $N=16$ modes, and with varying $\mu=0,10,100$ (we also ran $\mu=1$, and results are very close to those for $\mu=0$). Results are shown in figure 3 for energy and drag prediction, and we observe a big improvement from using DA. For $N=16$ and $\mu=100$, very good accuracy is achieved from the DA-ROM. For $N=8$, $\mu=10$ is somewhat more accurate than for $\mu=100$, but both are better than no DA.

4.3 Inaccurate Snapshots

In this section, we investigate the DA-ROM performance when the snapshots are inaccurate. Specifically, we consider the same test as in Section 4.2, but now with only a small amount of data being used to build the ROM basis. This is an important aspect of the ROM to investigate, because in practical applications complete data is generally not available, or the amount of data needed to sufficiently capture the behavior of the true solution is unknown.

We generated these inaccurate snapshots for $Re$=500 using less than one period of data: basis 1 used 64% of one period of data while basis 2 used 84%. See figure 4 for the first five basis functions generated by the ROM; the basis functions for the full ROM are also included for comparison.

In figure 5, we show the results of the DA-ROM using only 8 modes, with basis 1 and 2 defined above, and $\mu$ ranging from 100 to 500. DA significantly improves the accuracy of the ROM, and basis 2 does better at predicting the drag coefficient than basis 1.

Figure 6 shows energy and drag coefficient plots versus time using $N=12$ modes, and the nudging parameter $\mu$ is varied from 100 to 500. We see similar results as the case of using 8 modes; for both bases, DA significantly improves the accuracy of the ROM, compared to the ROM without DA ($\mu=0$), which becomes more and more inaccurate as time goes on. Basis 2 is able to very accurately predict the drag coefficient.

The results in this section suggest that DA can dramatically improve the accuracy of a ROM when insufficient data is available to build the ROM, which is the general case in practical applications. We also emphasize that the improvement in the DA-ROM accuracy over the standard ROM accuracy is significantly larger in the realistic case of inaccurate snapshot construction. Indeed, comparing figures 5 and 6 with figure 3, we notice that the absolute improvement in the DA-ROM is much larger in the former than in the latter (this could be clearly seen from the magnitude of the $y$-axis).

4.4 Adaptive Nudging

To further improve the accuracy of the DA-ROM solution, we also consider nudging that is adaptive in time. While the error estimate we prove guarantees convergence up to discretization error and ROM truncation error exponentially fast in time, it may not be sufficient to expect good numerical results. In practice, the ROM truncation error is often quite large, and can make the error bounds be too large to guarantee accurate predictions, especially over long time intervals. We propose below an adaptive nudging technique that will help produce better results by forcing the DA-ROM predicted energy to be more accurate.

4.4.1 Algorithm

In this section, we propose to change $\mu$ adaptively in time, based on the accuracy of the energy prediction of the ROM as well as the sign of the contribution of the data assimilation term to the energy balance. The semi-discrete algorithm reads: Find $u_{r}\in X_{r}$ such that for all $v_{r}\in X_{r}$,

[TABLE]

with $v_{0}=P_{r}(u_{0})$, and $\mu$ is the adaptive nudging parameter.

We begin the discussion with an energy estimate. Choosing $\chi_{r}=u_{r}$ vanishes the nonlinear term, and after bounding the forcing term in the usual way we obtain the energy estimate

[TABLE]

We assume this estimate is sharp in the following analysis, and that we know $\|u(t^{n})\|$ in addition to $I_{H}(u)(t^{n})$.

The adaptive strategy is to adjust $\mu$ so the contribution of the data assimilation term removes dissipation if the ROM-DA energy is too small, and adds dissipation if the energy is too large. We use the term dissipation loosely, since here we refer to dissipation from the DA term only meaning that it adds positivity to the left hand side of the energy estimate. Now after step $n$ we can calculate (1) the DA-ROM energy $\frac{1}{2}\|u_{r}^{n}\|^{2}$ and the true energy $\frac{1}{2}\|u(t^{n})\|^{2}$; and (2) the sign of the contribution of the data assimilation term (DAT):

[TABLE]

With this information, we check the energy error to see if it is too high (or too low), and if so, then add dissipation by increasing $\mu$ if $DAT>0$ and decreasing $\mu$ otherwise; or do the opposite to decrease dissipation.

How often to adjust $\mu$, and by how much each time, are interesting questions. In our numerical tests below, we checked the value of $DAT$ every 10 time steps, since there is some calculation cost involved, and changed $\mu$ by $\pm 1$ each time, as large sudden changes in $\mu$ gave bad results.

4.4.2 Numerical Results

We follow the same problem set up outlined in Section 4.2 (again using the full ROM basis), but now choosing $\mu$ adaptively in time. We note that, in addition to the Reynolds number we considered in the previous numerical experiments (i.e., $Re=500$), we also consider $Re=1000$. Figures 7 and 8 show the energy and drag plots for the DA-ROM algorithm with the adaptive nudging described above, and for constant $\mu$, for no DA. For both $Re$, the adaptive DA-ROM yields the most accurate results, outperforming the ROM without DA, and the DA-ROM with a constant $\mu$. Also included are plots of the $\mu$ values chosen by the algorithm at each timestep. We observe that the behavior of the values of $\mu$ is similar to that of the plots of $DAT$ in the figures.

5 Conclusions

In this paper, we put forth a new data assimilation reduced order model (DA-ROM) for fluid flows. The new DA-ROM adds to the standard ROM a feedback control term that nudges the ROM approximation towards the reference solution corresponding to the observed data. The new DA-ROM’s implementation is extremely simple: The nudging term can be implemented into existing codes completely at the linear algebraic level, without any changes to the rest of the discretization. The nudging term dramatically increases the accuracy of the new DA-ROM by utilizing the available low-resolution data, without the need to increase the number of ROM basis functions. We proved that with a properly chosen nudging parameter, the new DA-ROM algorithm converges exponentially fast in time to the true solution, up to discretization and ROM truncation errors. We also proposed a strategy for nudging adaptively in time, by adding or removing dissipation arising from the nudging to better match true solution energy. Finally, we performed a numerical investigation of the new DA-ROM in the simulation of a 2D flow past a circular cylinder. The numerical results showed that the adaptive nudging DA-ROM significantly improves the long time ROM accuracy, especially when the snapshots used to construct the ROM are inaccurate, which is generally the case in realistic applications.

We intend to pursue several research avenues. First, we will investigate whether numerical analysis can help determine the optimal parameter in the adaptive nudging approach for the new DA-ROM. We also want to extend the numerical investigation of the DA-ROM to complex 3D flows. Finally, we will examine whether using a spectral type of nudging in the DA-ROM instead of the current physical nudging yields better results.