A non-autonomous time-delayed SIR model for COVID-19 epidemics prediction in China during the transmission of Omicron variant

1.   Introduction

Performing accurate simulations of complex fluid flows that match real-world observations or experiments typically requires highly precise knowledge of the initial data. However, such data is often known in very sparsely-distributed locations, which is the case in, e.g., weather observation, ocean monitoring, etc. Thus, accurate, deterministic simulations based on initial data are often impractical. Data assimilation is a collection of methods that works around this difficulty by incorporating incoming data into the simulation to increase accuracy, hence data assimilation techniques are highly desirable to incorporate into simulations. However, the underlying physical equations often suffer from stability issues which can reduce the accuracy gained by using data assimilation. While there are many ways to stabilize numerical simulations, it is far from obvious how to adapt data assimilation techniques to combine them with cutting-edge stabilization methods. Therefore it becomes worthwhile to seek new ways to incorporate data assimilation into stabilized schemes. In this article, we propose and analyze a new approach to this problem which combines continuous data assimilation with velocity-vorticity stabilization.

Since Kalman's seminal paper [29] in 1960, a wide variety of data assimilation algorithms have arisen (see, e.g., [14,32]). In [5], Azouani, Olson, and Titi proposed a new algorithm known as continuous data assimilation (CDA), also referred to as the AOT algorithm. Their approach revived the so-called "nudging" methods of the 1970's (see, e.g., [4,26]), but with the addition of a spatial interpolation operator. This seemingly minor change had profound impacts, and the authors of [5] were able to prove that using only sparse observations, the CDA algorithm applied to the 2D Navier-Stokes equations converges to the correct solution exponentially fast in time, independent of the choice initial data. This stimulated a large amount of recent research on the CDA algorithm; see, e.g., [3,6,7,10,11,13,17,18,19,20,21,22,27,31,30,35,40,41] and the references therein. The recent paper [15] showed that CDA can be effectively used for weather prediction, showing that it can indeed be a powerful tool on practical large scale problems. Convergence of discretizations of CDA models was studied in [30,41,27,21], and found results similar to those at the continuous level. Our interest in the CDA algorithm arises from its adaptability to a wide range of nonlinear problems, as well as its small computational cost and straight-forward implementation. These qualities make it an ideal candidate for combining data assimilation with stabilization techniques; in particular, with the recently developed velocity-vorticity stabilization, described below.

Flows of incompressible, viscous Newtonian fluids are modeled by the Navier-Stokes equations (NSE), which take the form

together with suitable boundary and initial conditions. Here, $u$ denotes a velocity vector field, $p$ is pressure, $f$ is external (given) force, and $\nu>0$ represents the kinematic viscosity which is inversely proportional to the Reynolds number. Solving the NSE is important in many applications, however it is well known that doing so can be quite difficult, especially for small $\nu$. Many different tools have been used for more accurate numerical simulations of the NSE, for example using NSE formulations tailored to particular application problems [23,12,37,33] or discretization and stabilization methods [39,28], and more recently using observed data to improve simulation [5,30,43,44,8].

We consider in this paper discretizations of a continuous data assimilation (CDA) enhancement applied to the following velocity-vorticity (VV) formulation of the 2D NSE:

Here, $\omega$ represents the (scalar) vorticity, $P: = p+\tfrac12|u|^2$ is the Bernoulli pressure, and rot is the 2D curl operation: ${\rm{rot}} \binom{f_1}{f_2}: = \tfrac{\partial f_2}{\partial x} - \tfrac{\partial f_1}{\partial y}$. In the NSE, the velocity and vorticity are coupled via the relationship $\omega = {\rm{rot}} u$ (or equivalently, the Biot-Savart Law). However, the VV formulation typically does not enforce this relationship, so $u$ and $\omega$ are only coupled via the evolution equations in (2), and the relationship $\omega = {\rm{rot}} u$ is recovered a posteriori, so that at the continuous level, (2) is formally equivalent to (1). However, in practice, discretizations of VV can behave quite differently from typical discretizations of NSE, providing better stability as well as accuracy (especially for vorticity) for vortex dominated or strongly rotating flows, see [36,37,34,2] and references therein. A very interesting property of (2) was recently shown in [25], where it was proven that the system (2) when discretized with standard finite elements and a decoupling backward Euler or BDF2 temporal discretization was unconditionally long-time stable in both $L^2$ and $H^1$ norms for both velocity and vorticity; no such analogous result is known for velocity-pressure discretizations/schemes. Hence the scheme itself is stabilizing, even though it is still formally consistent with the NSE. The recent work in [2] showed that these unconditionally long-time stable schemes also provide optimal vorticity accuracy, yielding a vorticity solution that is one full order of spatial accuracy better than for an analogous velocity-pressure scheme.

We consider herein CDA applied to (2), which yields a model of the form

where $I_H$ is an appropriate interpolation operator, $I_H(u)$ and $I_H(\omega)$ are assumed known from measurement data (i.e. $u$ and $\omega$ are known at some points in space), and $\mu_1,\mu_2\ge 0$ are nudging parameters. If $\mu_2 = 0$, then vorticity is not nudged and $I_H(\omega)$ need not be assumed known. Due to the success of (2) in recent papers [25,2,36] and that of CDA in the works mentioned above, combining these ideas and studying (2) is a natural next step to see whether CDA will provide optimal long-time accuracy for the VV schemes already known to be unconditionally long-time stable. Herein, we do find that CDA provides convergence of (3), with any initial condition, to the true NSE solution (up to optimal discretization error) and moreover that CDA preserves the long-time stability.

This paper is organized as follows. In Section 2, we introduce the necessary notation and preliminaries needed in the analysis. In Section 3, we propose and analyze a fully discrete scheme for (3), and show that for nudging velocity and vorticity together and nudging just velocity, algorithms are long-time stable in $L^2$ and $H^1$ norms and long-time optimally accurate in $L^2$ velocity and vorticity (under the usual CDA assumptions on the coarse mesh and nudging parameter). In Section 4, we illustrate the theory with numerical tests, and finally draw conclusions in section 5.

2.   Notation and preliminaries

We now provide notation and mathematical preliminaries to allow for a smooth analysis to follow. We consider the domain $\Omega \subset \mathbb{R}^2$ to be the $2\pi$-periodic box, with the $L^2(\Omega)$ norm and inner product denoted by $\| \cdot \|$ and $(\cdot, \cdot)$ respectively, while all other norms will be appropriately labeled.

For simplicity, we use herein periodic boundary conditions for velocity and vorticity. Extension to full nonhomogeneous Dirichlet conditions can be performed by following analysis in [34], although for no-slip velocity together with the more physically consistent natural vorticity boundary condition studied in [36,38] more work would be needed to handle the boundary integrals. We denote the natural corresponding function spaces for velocity, pressure, and vorticity by

In $X$ (and $W$), we have the Poincaré inequality: there exists a constant $C_P$ depending only on $\Omega$ such that for any $\phi\in X$ (or $W$),

We define the skew-symmetric trilinear operator $b^*:X\times W \times W \rightarrow \mathbb{R}$ to use for the nonlinear term in the vorticity equation, by

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

Lemma 2.1. Suppose constants $r$ and $B$ satisfy $r>1$, $B\ge 0$. Then if the sequence of real numbers $\{a_n\}$ satisfies

we have that

2.1.   Discretization preliminaries

Denote by $\tau_h$ a regular, conforming triangulation of the domain $\Omega$, and let $X_h \subset X$, $Q_h \subset Q$ be velocity-pressure spaces that satisfy the inf-sup condition. We will assume the use of $X_h = X \cap P_k(\tau_h)$ and $Q_h = Q \cap P_{k-1}(\tau_h)$ Taylor-Hood or Scott-Vogelius elements (on appropriate meshes and/or polynomial degrees, see [24] and references therein). The discrete vorticity space is defined as $W_h : = W \cap P_k(\tau_h).$ Define the discretely divergence free subspace by

We will assume the mesh is sufficiently regular so that the inverse inequality holds in $X_h$: There exists a constant $C$ such that

The discrete Laplacian operator is defined as: For $\phi\in H^1(\Omega)^2$, $\Delta_h \phi \in X_h$ satisfies

The definition for $\Delta_h$ is written the same way when applied in $W_h$, since this is simply the above definition restricted to a single component.

The discrete Stokes operator $A_h$ is defined as: For $\phi\in H^1(\Omega)^2$, find $A_h \phi \in V_h$ such that for all $v_h \in V_h$,

By the definition of discrete Laplace and Stokes operators, we have the Poincaré inequalities

We recall the following discrete Agmon inequalities and discrete $L^p$ bounds [25]:

We note that all bounds above for $X_h$ trivially hold in $W_h$, since $W_h$ functions can be considered as components of functions in $X_h$.

A function space for measurement data interpolation is also needed. Hence we require another regular conforming mesh $\tau_H$, and define $X_H = P_r(\tau_H)^2$ and $W_H = P_r(\tau_H)$ for some polynomial degree $r$. We require that the coarse mesh interpolation operator $I_H$ used for data assimilation satisfies the following bounds: for any $w \in H^1(\Omega)^d$,

These are key properties for the interpolation operator that allow for both mathematical theory as well as providing guidance on how small $H$ should be (i.e. how many measurement points are needed). We note the same $I_H$ operator is used for vector functions and scalar functions, with it being applied component-wise for vector functions.

3.   Analysis of a CDA-VV scheme

We consider now a discretization of (3) that uses a finite element spatial discretization and backward Euler temporal discretization. The backward Euler discretization is chosen only for simplicity of analysis; all results extend to the analogous BDF2 scheme following analysis in [2,25]. One difference of our scheme below compared to other discretizations of CDA is that $I_H$ is also applied to the test functions in the nudging terms. This was first proposed by the authors in [41], and allows for a simpler stability analysis as well as to the use of special types of efficient interpolation operators.

Algorithm 3.1. Given $v_h^0\in V_h$ and $w_h^0\in W_h$, find $(v_h^{n+1}, w_h^{n+1}, P_h^{n+1})$ $\in$ $(X_h, W_h, Q_h)$ for $n = 0,1,2,...$, satisfying

for all $(\chi_h, \psi_h, r_h) \in (X_h,W_h,Q_h)$, where $I_H(u^{n+1})$, $I_H({\rm{rot}} u^{n+1})$ are assumed known for all $n\ge 0$.

We begin our analysis with long-time stability estimates, followed by long-time accuracy.

3.1.   Stability analysis of Algorithm 3.1

In this subsection, we prove that Algorithm 3.1 is unconditionally long-time $L^2$ and $H^1$ stable for both velocity and vorticity. This property was proven for the scheme without nudging in [25], and so these results show that CDA preserves this important property that is (seemingly) unique to VV schemes of this form.

Lemma 3.2 ($L^2$ stability of velocity and vorticity). Let $f \in L^\infty(0,\infty; L^2)$ and $u\in L^{\infty}(0,\infty; H^1)$. Then, for any ${\Delta t}>0$, any integer $n>0$, and nudging parameters $\mu_1,\mu_2\ge 0$, velocity and vorticity solutions to Algorithm 3.1 satisfy

where $\alpha = 1+\nu C_P^{-2}\Delta t$.

Proof. Begin by choosing $\chi_h = 2\Delta t v_h^{n+1}$ in (13), which vanishes the nonlinear and pressure terms, and leaves

The first right hand side term is bounded using the dual norm and Young's inequality via

and for the interpolation term, we use Cauchy-Schwarz, the interpolation property (12) and Young's inequality to get

Combining the above estimates and dropping $\|v_h^{n+1}-v_h^{n}\|^2$ and $\|I_H(v^{n+1}_{h})\|^2$ from the left hand side produces the bound

and thanks to the Poincaré inequality, we obtain

Defining $\alpha = 1+ C_P^{-2} \Delta t \nu> 1$, and then applying Lemma 2.1 reveals the $L^2$ stability bound (16) for the velocity solution of Algorithm 3.1.

Applying similar analysis to the above will produce the stated $L^2$ vorticity bound.

Lemma 3.3 ($H^1$ stability of velocity and vorticity). Let $f \in L^\infty(0,\infty; H^1)$ and $u \in L^{\infty}(0,\infty; H^1)$. Then, for any ${\Delta t}>0$, any integer $n>0$, and nudging parameters $\mu_1,\mu_2\ge 0$, velocity and vorticity solutions to Algorithm 3.1 satisfy

where $\alpha = 1+\nu C_P^{-2}\Delta t$.

Proof. After testing the velocity equation (13) with $\chi_h = 2 {\Delta t} A_hv_h^{n+1}$, we obtain

We now bound the right hand side terms. First, the forcing term is bounded by Cauchy-Schwarz and Young's inequalities via

Then, for the nonlinear terms, we again apply Hölder, discrete Agmon (9) and generalized Young inequalities, and the result of Lemma 3.2 to get

Lastly, the interpolation term is bounded using Cauchy-Schwarz and interpolation property 12, followed by Young's inequality and the result of Lemma 3.2 to obtain

Combining all these bounds for right hand side terms and dropping nonnegative term $\|\nabla(v_h^{n+1}-v_h^{n})\|^2$ on left hand side give us that

By the Poincaré inequality (7), we now get

where $\alpha = 1+\nu C_P^{-2} \Delta t$. Finally, we apply Lemma 2.1 and reveal (18).

For the vorticity estimate, choose $\psi_h = 2 {\Delta t} \Delta_h w_h^{n+1}$ in (15) to get

From here, the proof follows the same strategy as the $H^1$ velocity proof above, except the nonlinear term is handled slightly differently. We use the discrete Agmon inequality (9), the discrete Sobolev inequality (10), the result of Lemma 3.1, the $H^1$ stability bound for vorticity (18) proven above, and the generalized Young's inequality, as follows.

Now proceeding as in the velocity $H^1$ bound will produce the $H^1$ vorticity stability bound (19).

3.2.   Long-time accuracy of Algorithm 3.1

We now consider the difference between the solutions of (13) - (15) to the NSE solution. We will show that the algorithm solution converges to the true solution, up to an optimal $O( {\Delta t} + h^{k+1})$ discretization error, independent of the initial condition, provided a restriction on the coarse mesh width and nudging parameters. We will give two results, the first for $\mu_2>0$ and the second for $\mu_2 = 0$; while they both provide optimal long-time accuracy, when $\mu_2>0$ the convergence to the true solution occurs more rapidly in time.

In our theory below for long-time accuracy of Algorithm 3.1, we assume the use of Scott-Vogelius elements. This is done for simplicity, as for non-divergence-free elements like Taylor-Hood elements, similar optimal results can be obtained (although with some additional terms and different constants) but require more technical details; see, e.g., [21]. We define the following projection operator, $P_V$, which will be used in the following accuracy analysis: Given $\phi \in H^1(\Omega)$ with $\nabla \cdot \phi = 0$, $P_V(\phi)\in V_h$ satisfies $(\nabla (\phi - P_V(\phi)),\nabla v_h) = 0$ for all $v_h\in V_h$. We also define $P_W: H^1(\Omega) \rightarrow W_h$ in the same way, with $V_h$ replaced by $W_h$. The operators $P_V$ and $P_W$ are known to have optimal approximation properties on the divergence free subspace of $X$ and on $W$, respectively, in both $L^2$ and $H^1$ norms, provided some commonly assumed properties of the domain [9,42].

Theorem 3.4 (Long-time $L^2$ accuracy of Algorithm 3.1 with $\mu_1>0$ and $\mu_2>0$). Let true solutions $u \in L^\infty(0, \infty ; H^{k+2}(\Omega))$, $p \in L^\infty(0, \infty ; H^{k}(\Omega))$ where $k\geq 1$ and $u_{t}, u_{tt}\in L^\infty (0, \infty ; H^1)$. Then, assume that time step ${\Delta t}$ is sufficiently small, and that $\mu_1$ and $\mu_2$ satisfy

where $H$ is chosen so that this inequality holds. Then, for any time $t^n$, $n = 0, 1, 2,...$, we have for solutions of Algorithm 3.1 using Scott-Vogelius elements,

where

and $\lambda = \min \left\{ \frac{\mu_1}{4} + \frac{\nu C_P^{-2}}{4} , \frac{\mu_2}{4} + \frac{\nu C_P^{-2}}{4} \right\}$ with $C$ independent of ${\Delta t}$, $h$ and $H$.

Proof. The true NSE solution satisfies the VV system

where $u^n$ is the velocity at time $t^n$, $P^n$ the Bernoulli pressure, $\omega^n : = {\rm{rot}} u^n$, and $t^*, t^{**} \in [t^n, t^{n+1}]$. Note that by Taylor expansion, we can write $\omega^{n}-\omega^{n+1} = - {\Delta t} \omega_t (s^*)$ where $s^* \in [t^n, t^{n+1}]$.

The difference equations are obtained by subtracting the solutions to Algorithm 3.1 from (22)-(24) after testing them with test functions from $\chi_h\in X_h$, $r_h\in Q_h$ and $\psi_h\in W_h$, respectively. We define the differences between velocity and vorticity as $e_v^{n} : = v_h^n - u^n$ and $e_w^n : = w_h^n - \omega^n$, respectively. Next, we will decompose the error into a term that lies in the discrete space $V_h$ and one outside. To do so, add and subtract the discrete Stokes projection of $u^n$, denoted $P_V u^{n}$, to $e_v^n$ and let $\eta_v^n : = P_V u^{n} - u^n$, $\phi_{h,v}^n : = v_h^n - P_V u^{n}$. Then $e_v^n = \phi_{h,v}^n + \eta_v^n$ and $\phi_{h,v}^{n} \in V_h$. In a similar manner, by taking the $H^1_0$ projection of ${\rm{rot}} u^n$ into $W_h$ denoted by $P_W \omega^{n}$, we obtain $e_w^n = \phi_{h,w}^n + \eta_w^n$ with $\phi_{h,w}^n \in V_h$.

For velocity, since $(\nabla \eta_v^{n+1},\nabla \phi_{h,v}^{n+1}) = 0$, the difference equation becomes

and similarly for vorticity, we have

where in (25) we have added and subtracted $\phi_{h,v}^{n+1}$ to write it in the form found above using

and similarly for (26).

Next, we bound the terms on right hand side of difference equations, starting with the velocity difference equation (25). The first three right hand side terms are bounded using Cauchy-Schwarz and Young's inequalities, via

where $s^* \in [t^n, t^{n+1}]$.

For nonlinear terms in (25), first we add and subtract $e_w^{n+1}$ in the first component, and $u^{n+1}$ in second component to get

The all resulting terms are bounded by Hölder's and Young's inequalities to obtain

Then, for the last nonlinear term, we apply Hölder's, Poincaré's and Young's inequalities and get

Next, the first interpolation term on the right hand side of (25) will be bounded with Cauchy-Schwarz inequality and (11) to obtain

For the second interpolation term, we apply inequality (12), which yields

Finally, the last interpolation term will be bounded using Cauchy-Schwarz, (12), and Young's inequality to get the bound

We now move on to the vorticity difference equation, (26). All the linear terms are majorized in a similar manner as in the velocity case, and so we show below the bounds only for the nonlinear terms. Due to the use of Scott-Vogelius elements, the skew-symmetric form reduces to the usual convective form, so $b^*(u,v,w) = (u\cdot\nabla v,w)$, with $\| \nabla \cdot u \| = 0$. To bound the first nonlinear term on the right hand side of (26), we begin by breaking up the velocity error term, then apply Hölder's and Young's inequalities, yielding

For the second nonlinear term, we use Hölder's, Póincare's and Young's inequalities, which gives

For the last nonlinear term, we begin by breaking up the velocity error term, then apply Hölder's and Young's inequalities to get

Replacing the right hand sides of (25) and (26) with the computed bounds and dropping nonnegative terms with $\| \phi_{h,v}^{n+1}- \phi_{h,v}^{n}\|^2$ yields the bound

Using the assumptions on $H$ and the nudging parameters, the time step restriction, and smoothness of the true solution, this reduces to

Now define

Using this in the inequality after applying Poincare's inequality and multiplying each side by $2 {\Delta t}$, we get

Then, with $R: = \left(\mu_1^{-1} {\Delta t}^2 + \nu^{-1} {\Delta t}^2 + \nu^{-1} h^{2k+2} + \mu_1 h^{2k+2} + \mu_2 h^{2k+2}\right)$ and $\lambda : = \min \left\{ \lambda_1, \lambda_2 \right\}$, we obtain the bound

By Lemma 2.1, this implies

Lastly, applying triangle inequality completes the proof.

Theorem 3.5 (Long-time $L^2$ accuracy of Algorithm 3.1 with $\mu_1>0, \mu_2 = 0$). Let true solution $u \in L^\infty(0, \infty ; H^{k+2}(\Omega))$, $p \in L^\infty(0, \infty ; H^{k}(\Omega))$ where $k\geq 1$ and $u_{t}, u_{tt}, \in L^\infty (0, \infty ; H^1)$. Then, assume that time step ${\Delta t}$ is sufficiently small, $\mu_2 = 0$, and that $\mu_1$ satisfies

where $H$ is chosen so that this inequality holds. Then for any time $t^n$, $n = 0, 1, 2,...$, solutions of of Algorithm 3.1 using Scott-Vogelius elements satisfy

where

and $\lambda = \min \left\{ \frac{\mu_1}{4}+\frac{C_P^{-2} \nu}{4}, \frac{C_P^{-2}\nu}{4}\right\}$ with $C$ independent of ${\Delta t}$, $h$ and $H$.

Remark 1. Algorithm 3.1 converges to the true solutions up to optimal discretization error in both cases $\mu_1, \mu_2 >0$ and $\mu_1>0, \mu_2 = 0$. The key difference between two cases is that when $\mu_2 = 0$, the convergence in time to reach optimal accuracy is much slower since $\lambda$ does not scale with the nudging parameters. This phenomena is illustrated in our numerical tests.

Proof. We follow the same steps with the proof of Theorem 3.4. The difference equation for velocity is already the same with (25), and just two nonlinear terms in the velocity difference equation are bounded with differently in this case. By Hölder, Poincaré and Young's inequalities, we get the bounds

All terms on the right hand side of vorticity difference equation for Theorem 3.4 are bounded identically. Proceeding as in the previous proof, we arrive at

Provided ${\Delta t}$ is sufficiently small and the restriction

holds, then applying Poincaré inequality to the terms on left hand side and using

and assumptions on the true solution, we obtain

From here, the proof is finished in the same way as the previous theorem.

3.3.   Second order temporal discretization

We now present results for a second order analogue of the first order algorithm studied above.

Algorithm 3.6. Find $(v_h^{n+1}, w_h^{n+1}, q_h^{n+1}) \in (X_h, W_h, Q_h)$ for $n = 0,1,2,...$, satisfying

for all $(\chi_h, \psi_h, r_h) \in (X_h,W_h,Q_h)$, with $v^0\in X$ and $I_H(u^{n+1})$, $I_H({\rm{rot}} u^{n+1})$ given.

Stability and convergence results follow in the same manner as the first order scheme results above, using G-stability theory as in [1,2,30].

Theorem 3.7 (Long-time stability and accuracy of Algorithm 3.6 with $\mu_1>0$ and $\mu_2>0$). For any time step ${\Delta t}>0$, and any time $t^n$, $n = 0, 1, 2,...$, we have that solutions of Algorithm 3.6 satisfy

with $C$ independent of $n$, ${\Delta t},\ h,\ H$.

Furthermore, if we suppose the true solution $u \in L^\infty(0, \infty ; H^{k+2}(\Omega))$, $p \in L^\infty(0, \infty ; H^{k}(\Omega))$ where $k\geq 1$ and $u_{t}, u_{ttt}, \in L^\infty (0, \infty ; H^1)$, that time step ${\Delta t}$ is sufficiently small, Scott-Vogelius elements are used, and that $\mu_1$ and $\mu_2$ satisfy $C(u) \leq \mu_1,\mu_2 \leq \frac{C\nu}{H^2}$, we have the bound

where

and $\lambda = \min \left\{ \frac{\mu_1}{4} + \frac{\nu C_P^{-2}}{4} , \frac{\mu_2}{4} + \frac{\nu C_P^{-2}}{4} \right\}$ with $C$ independent of ${\Delta t}$, $h$ and $H$.

4.   Numerical experiments

In this section, we illustrate the above theory with two numerical tests, both using Algorithm 3.6. Our first test is for convergence rates on a problem with analytical solution, and the second test is for flow past a flat plate. For both tests, we report results only for $(P_2,P_1)$ Taylor-Hood elements for velocity and pressure, and $P_2$ for vorticity; however we also tried Scott-Vogelius elements on barycenter refined meshes that produced similar numbers of degrees of freedom, and results were very similar to those of Taylor-Hood. The coarse velocity and vorticity spaces $X_H$ and $W_H$ are defined to be piecewise constants on the same mesh used for the computations. The interpolation operator $I_H$ was taken to be the $L^2$ projection operator onto $X_H$ (or $W_H$), which is known to satisfy (11)-(12) [16].

4.1.   Experiment 1: Convergence rate test

For our first test, we investigate the theory above for Algorithm 3.6. Here we use the analytic solution

on the unit square domain $\Omega = (0,1)^2$ with kinematic viscosity $\nu = 1.0$, and use the NSE to determine $f$ and boundary conditions. We take the final time $T = 1$, and choose initials conditions for Algorithm 3.6's velocity and vorticity to be 0. For the discretization, $(P_2,P_1)$ Taylor-Hood elements are used for velocity and pressure, $P_2$ for vorticity, and a time step size of $\Delta t = 0.001$. From Section 3, we expect third order spatial convergence rate in the $L^2$ norm for large enough times. Results are presented below for two cases, $\mu_2>0$ and $\mu_2 = 0$.

4.1.1.   Results for $\mu_1>0$ and $\mu_2>0$

To test this case, we first calculated spatial convergence rates at the final time $T = 1$ with the $L^2$ error, using successively refined uniform meshes and $\mu_1 = \mu_2 = 100$. Errors and rates are shown in table 1, and show clear third order spatial convergence of both velocity and vorticity. Deterioration of the rates for the smallest $h$ is expected since the time step $\Delta t$ is fixed while the spatial mesh width decreases.

We next consider convergence to the true solution exponentially in time (up to discretization error). Here we take $h = 1/32$, and compute solutions using $\mu_1 = \mu_2 = \mu$, with $\mu = 1,\ 10,\ 100,\ 1000$. Results are shown in figure 1, as $L^2$ error versus time for velocity and vorticity. We observe exponential convergence in time of both velocity and vorticity, up to about $10^{-5}$, which is consistent with the choices of $h$ and $\Delta t$ and the accuracy of the method. We note that as $\mu$ is increased, convergence is faster in time, which is consistent with our theory for the case of $\mu_1>0$ and $\mu_2>0$.

4.1.2.   Results for $\mu_1>0$ and $\mu_2=0$

We now consider the same tests as above, but with $\mu_2 = 0$. This is an important case, since it is not always practical to obtain accurate vorticity measurement data. Just as in the first case, we first calculated spatial convergence rates at the final time $T = 1$ for the $L^2$ error, on the same successively refined uniform meshes, but now with $\mu_1 = 100$ and $\mu_2 = 0$. Errors and rates are shown in table 2, and show clear third order spatial convergence of both velocity and vorticity. Deterioration of the rates for the smallest $h$ is expected since the time step $\Delta t$ was fixed at 0.001, although the vorticity errors are slightly worse than for the case of $\mu_2 = 100$ shown in table 1, and the deterioration of the rates occurs a bit earlier. Hence we observe essentially the same velocity errors and rates compared to the case of $\mu_2 = 100$, and slightly worse vorticity error but still with optimal $L^2$ accuracy.

To test exponential convergence in time for the case of $\mu_2 = 0$, we again take $h = 1/32$, and compute solutions using $\mu_1 = 1,\ 10,\ 100,\ 1000$. Results are shown in figure 2, as $L^2$ error versus time for velocity and vorticity. Although we do observe exponential convergence in time of both velocity and vorticity, up to about $10^{-5}$ which is the same accuracy reached when $\mu_2 = 100$ above. An important difference here compared to when $\mu_2 = 100$ is that the convergence of vorticity to the true solution is independent of $\mu_1$, and the convergence of velocity is slower for larger choices of $\mu_1$. This reduced dependence of the convergence on the nudging parameters when $\mu_2 = 0$ is consistent with our theory. Hence without vorticity nudging, long-time optimal accuracy is still achieved, but it takes longer in time to get there.

4.2.   Experiment 2: Flow past a normal flat plate

To test Algorithm 3.6 on a more practical problem, we consider flow past a flat plate with $Re = 50$. The domain of this problem is $[-7,20]\times[-10,10]$ rectangular channel with a $0.125\times1$ plate fixed ten units into the channel from the left, vertically centered. The inflow velocity is $u_{in} = \langle0,1\rangle$, no-slip velocity and the corresponding natural vorticity boundary condition from [36] are used on the walls and plate, and homogeneous Neumann conditions are enforced weakly at the outflow. A setup for the domain is shown in figure 3. There is no external forcing applied, $f = 0$. The viscosity is taken to be $\nu = 1/50$ which is inversely proportional to $Re$, based on the height of the plate. The end time for the test is $T = 80$. A DNS was run until for 160 time units (from t = -80 to t = 80), and for $t>0$ measurement data for the VV-DA simulation was sampled from the DNS.

We computed solutions using a Delaunay generated triangular meshes that provided $27,373$ total degree of freedom with $(P_2, P_1, P_2)$ velocity-pressure-vorticity elements, and time step ${\Delta t} = 0.02$. We first compared convergence in time to the DNS solution, for two cases: $\mu_1 = \mu_2>0$ and $\mu_1>0,\ \mu_2 = 0$. Plots of $L^2$ velocity and vorticity error for both of these cases are shown in figure 4, with varying nudging parameters. There is a clear advantage seen in the plots for the simulations with $\mu_2>0$: when vorticity is nudged in addition to velocity, convergence to the true solution is much faster in time. The convergence when $\mu_2 = 0$ appears to still be occurring, but is much slower and even by $t = 80$ the $L^2$ vorticity error is barely smaller than $O(1)$. We note that just like in the analytical test problem, when $\mu_2 = 0$ the vorticity convergence in time is independent of $\mu_1$.

To further illustrate the convergence of the DNS, we show contour plots of the DNS solution, the VVDA solution, and their difference, in figures 5-8. For these simulations, we used $\mu_1 = \mu_2 = 10$ in figures 5-6, but used $\mu_2 = 0$ for figures 7-8. As expected due to the plots in figure 4, when $\mu_1 = \mu_2 = 10$ we observe rapid convergence in the plots for VV-DA velocity and vorticity to the DNS velocity and vorticity, and by $t = 1$ the contour plots are visually indistinguishable. This is not the case, however, when $\mu_2 = 0$. In this case, while the velocity plots do agree with DNS velocities by $t = 1$ (in the eyeball norm), the vorticity error remains observable at $t = 10$ and even at $t = 20$ there are some small difference. The contour plots of the errors at early times for vorticity show the largest errors occur near vortex centers, indicating that the VV-DA method is not accurately predicting the strength of the vortices.

4.2.1.   Re=100

We also tested $Re = 100$ for flow past a flat plate, using the same discretization parameters as above for the $Re = 50$ case, and overall see similar results as for the $Re = 50$ case. When both velocity and vorticity are nudged, convergence to the true solution is exponential for both velocity and vorticity in the $L^2(\Omega)$ norm, and we observe that at early times the larger $\mu$ convergence curves are steeper, but at later times $\mu = 1$ converges more rapidly. When only velocity is nudged, larger $\mu$ provides faster convergence of the velocity, but the vorticity converges nearly independent of $\mu$.

5.   Conclusions

We have analyzed a VV scheme for NSE enhanced with CDA, using linearized backward Euler or BDF2 in time and finite elements in space. We proved that applying CDA preserves the unconditional stability properties of the scheme, and also yields optimal long-time accuracy if both velocity and vorticity are nudged, or velocity-only. If only velocity is nudged, then the convergence in time to the true solution is slower, but still exponentially fast in time. Numerical tests illustrate the theory, including the difference between nudging only velocity or also nudging vorticity.

For future directions, since nudging vorticity is difficult in practice due to accurate measurement data not typically being available, one may try to obtain better results for the velocity-only-nudging by penalizing the difference between $w_h$ and ${\rm{rot}} v_h$ in the vorticity equation. That is, by setting $\mu_2 = 0$ and adding the term $\gamma (w - {\rm{rot}} u)$ to the vorticity equation (3), it may be possible to analytically prove a convergence result resembling Theorem 3.4. Determining whether this is possible, and if so then for what values of $\gamma$, and whether it works in practice (i.e. how large are associated constants), would need a detailed further study which the authors plan to undertake.