Doubly-Adaptive Artificial Compression Methods for Incompressible Flow

TL;DR

This paper introduces adaptive artificial compression algorithms for incompressible flow that independently adjust time-step and compression parameters, supported by analysis and numerical validation, with minimal added complexity.

Contribution

It presents novel adaptive algorithms that independently tune time-step and artificial compression parameters, supported by theoretical analysis and numerical experiments.

Findings

Adaptive methods effectively control errors in incompressible flow simulations.

The algorithms maintain low computational complexity comparable to simpler methods.

Both first and second-order schemes are successfully embedded and tested.

Abstract

This report presents adaptive artificial compression methods in which the time-step and artificial compression parameter $\varepsilon $ are independently adapted. The resulting algorithms are supported by analysis and numerical tests. The first and second-order methods are embedded. As a result, the computational, cognitive and space complexities of the adaptive $% \varepsilon ,k$ algorithms are negligibly greater than that of the simplest, first-order, constant $\varepsilon ,$ constant $k$ artificial compression method.

Equations

$$\{\begin{array}{l}{u}_t+u\cdot \nabla u+\frac{1}{2}(\nabla \cdot u)u+\nabla p-\nu \Delta u=f\\ \epsilon p_t+\nabla \cdot u=0,\text{ where }0<\epsilon \text{ is small.}\end{array}$$

$$u^{\ast }=\left(1+\frac{k_{n+1}}{k_n}\right)u_n-\frac{k_{n+1}}{k_n}{u}_{n-1}\left(=2u_n-u_{n-1}\text{ for constant time-step}\right).$$

$$\frac{u_{n+1}-u_n}{k}+u^{\ast }\cdot \nabla u_{n+1}+\frac{1}{2}(\nabla \cdot u^{\ast })u_{n+1}+\nabla p_{n+1}-\nu \Delta u_{n+1}=f(t_{n+1}),$$

$$\epsilon \frac{p_{n+1}-p_n}{k}+\nabla \cdot u_{n+1}=0.$$

$$\frac{u_{n+1}-u_n}{k}+u^{\ast }\cdot \nabla u_{n+1}+\frac{1}{2}(\nabla \cdot u^{\ast })u_{n+1}-\frac{k}{\epsilon }\nabla \nabla \cdot u_{n+1}$$

$$-\nu \Delta u_{n+1}=-\nabla p_n+f_{n+1},$$

$$\text{ then given }{u}_{n+1}\text{: }{p}_{n+1}=p_n-\frac{k}{\epsilon }\nabla \cdot u_{n+1}.$$

$$\begin{array}[]{ccc}\text{GA-Method:}&&\frac{\varepsilon_{n+1}p_{n+1}-\sqrt{\varepsilon_{n+1}\varepsilon_{n}}p_{n}}{k_{n+1}}+\nabla\cdot u_{n+1}=0,\text{ or}\\ &&\\ \text{min-Method:}&&\frac{\varepsilon_{n+1}p_{n+1}-\min\{\varepsilon_{n+1},\varepsilon_{n}\}p_{n}}{k_{n+1}}+\nabla\cdot u_{n+1}=0.\end{array}$$

$$\frac{u_{n+1}-u_n}{k_{n+1}}+u^{\ast }\cdot \nabla u_{n+1}+\frac{1}{2}(\nabla \cdot u^{\ast })u_{n+1}+\nabla p_{n+1}-\nu \Delta u_{n+1}=f_{n+1}.$$

$$\begin{array}{rl}\frac{u_{n+1}^1-u_n}{k_{n+1}}& +u^{\ast }\cdot \nabla u_{n+1}^1+\frac{1}{2}(\nabla \cdot u^{\ast })u_{n+1}^1+\nabla p_{n+1}-\nu \Delta u_{n+1}^1=f_{n+1},\\ \text{For }\tau & =\frac{k_{n+1}}{k_n}\text{ let }{\alpha }_1=\frac{\tau (1+\tau )}{(1+2\tau )},\text{ then}:\\ u_{n+1}& =u_{n+1}^1-\frac{{\alpha }_1}{2}\left\{\frac{2k_n}{k_n+k_{n+1}}{u}_{n+1}^1-2u_n+\frac{2k_{n+1}}{k_n+k_{n+1}}{u}_{n-1}\right\}.\end{array}$$

$$D_2(n+1):=\frac{2k_n}{k_n+k_{n+1}}{u}_{n+1}^1-2u_n+\frac{2k_{n+1}}{k_n+k_{n+1}}{u}_{n-1}.$$

$$||v||={\left({\int }_{\Omega }|v(x){|}^{2}dx\right)}^{1/2}\text{ and }(v,w)={\int }_{\Omega }v(x)\cdot w(x)dx.$$

$$EST(1)=||u_{n+1}-u_{n+1}^1||=\frac{{\alpha }_1}{2}||D_2(n+1)||.$$

$$EST(2)$$

$${\alpha }_2$$

$$u_0=u_0(x)\text{ and }{p}_0=p_0(x),\text{ in }\Omega \text{,}$$

$$u_n=0\text{ on }\partial \Omega \text{ for }t>0\text{.}$$

$$\frac{u_{n+1}-u_n}{k_{n+1}}+u^{\ast }\cdot \nabla u_{n+1}+\frac{1}{2}(\nabla \cdot u^{\ast })u_{n+1}+\nabla p_{n+1}-\nu \Delta u_{n+1}=f_{n+1},$$

$$\frac{{\epsilon }_{n+1}{p}_{n+1}-\hat{\epsilon }{p}_n}{k_{n+1}}+\nabla \cdot u_{n+1}=0,\text{ where}$$

$$\hat{\epsilon }=\min \{{\epsilon }_{n+1},{\epsilon }_n\}\text{ for the min-Method and}$$

$$\hat{\epsilon }=\sqrt{{\epsilon }_{n+1},{\epsilon }_n}\text{ for the GA-Method }$$

$${\left({\epsilon }_{n+1}-{\epsilon }_n\right)}^{+}=\max \{0,{\epsilon }_{n+1}-{\epsilon }_n\}={\epsilon }_{n+1}-\min \{0,{\epsilon }_{n+1}-{\epsilon }_n\}.$$

$$\frac{1}{2}{\int }_{\Omega }|u_N{|}^2+{\epsilon }_N|p_N{|}^{2}dx+$$

$$\sum _{n=0}^{N-1}\frac{1}{2}{\int }_{\Omega }\min \{{\epsilon }_{n+1},{\epsilon }_n\}(p_{n+1}-p_n{)}^2+{\left({\epsilon }_{n+1}-{\epsilon }_n\right)}^{+}{p}_{n+1}{}^2$$

$$+{\left({\epsilon }_n-{\epsilon }_{n+1}\right)}^{+}{p}_n{}^{2}dx+\sum _{n=0}^{N-1}{\int }_{\Omega }\frac{1}{2}|u_{n+1}-u_n{|}^2+k_{n+1}\nu |\nabla u_{n+1}{|}^{2}dx$$

$$=\frac{1}{2}{\int }_{\Omega }|u_0{|}^2+{\epsilon }_0{p}_0{}^{2}dx+\sum _{n=0}^{N-1}{k}_{n+1}{\int }_{\Omega }{u}_{n+1}\cdot f_{n+1}dx.$$

$$\frac{1}{2}{\int }_{\Omega }|u_N{|}^2+{\epsilon }_N{p}_N{}^{2}dx+$$

$$\sum _{n=0}^{N-1}\frac{1}{2}{\int }_{\Omega }\min \{{\epsilon }_{n+1},{\epsilon }_n\}(p_{n+1}-p_n{)}^{2}dx+{\left({\epsilon }_{n+1}-{\epsilon }_n\right)}^{+}{p}_{n+1}^2$$

$$+{\left({\epsilon }_n-{\epsilon }_{n+1}\right)}^{+}{p}_n{}^{2}dx+\sum _{n=0}^{N-1}\frac{1}{2}{\int }_{\Omega }|u_{n+1}-u_n{|}^2+k_{n+1}\nu |\nabla u_{n+1}{|}^{2}dx$$

$$\le \frac{1}{2}{\int }_{\Omega }|u_0{|}^2+{\epsilon }_0{p}_0{}^{2}dx+\sum _{n=0}^{N-1}{k}_{n+1}\frac{1}{2\nu }||f_{n+1}|{|}_{-1}^2.$$

$$({\epsilon }_{n+1}{p}_{n+1}-\min \{{\epsilon }_{n+1},{\epsilon }_n\}{p}_n,p_{n+1})$$

$$={\epsilon }_{n+1}||p_{n+1}|{|}^2-\min \{{\epsilon }_{n+1},{\epsilon }_n\}(p_n,p_{n+1})$$

$$={\epsilon }_{n+1}||p_{n+1}|{|}^2-\min \{{\epsilon }_{n+1},{\epsilon }_n\}\left\{\frac{1}{2}||p_n|{|}^2+\frac{1}{2}||p_{n+1}|{|}^2-\frac{1}{2}||p_n-p_{n+1}|{|}^2\right\}$$

$$=\left({\epsilon }_{n+1}-\frac{1}{2}\min \{{\epsilon }_{n+1},{\epsilon }_n\}\right)||p_{n+1}|{|}^2$$

$$-\frac{1}{2}\min \{{\epsilon }_{n+1},{\epsilon }_n\}||p_n|{|}^2+\frac{1}{2}\min \{{\epsilon }_{n+1},{\epsilon }_n\}||p_n-p_{n+1}|{|}^2$$

$$=\frac{1}{2}{\epsilon }_{n+1}||p_{n+1}|{|}^2-\frac{1}{2}{\epsilon }_n||p_n|{|}^2+\frac{1}{2}\min \{{\epsilon }_{n+1},{\epsilon }_n\}||p_n-p_{n+1}|{|}^2+$$

$$+\frac{1}{2}\left({\epsilon }_{n+1}-\min \{{\epsilon }_{n+1},{\epsilon }_n\}\right)||p_{n+1}|{|}^2+\frac{1}{2}\left({\epsilon }_n-\min \{{\epsilon }_{n+1},{\epsilon }_n\}\right)||p_n|{|}^2.$$

$$({\epsilon }_{n+1}{p}_{n+1}-\min \{{\epsilon }_{n+1},{\epsilon }_n\}{p}_n,p_{n+1})=$$

$$=\frac{1}{2}{\epsilon }_{n+1}||p_{n+1}|{|}^2-\frac{1}{2}{\epsilon }_n||p_n|{|}^2+\frac{1}{2}\min \{{\epsilon }_{n+1},{\epsilon }_n\}||p_n-p_{n+1}|{|}^2+$$

$$+\frac{1}{2}{\left({\epsilon }_{n+1}-{\epsilon }_n\right)}^{+}||p_{n+1}|{|}^2+\frac{1}{2}{\left({\epsilon }_n-{\epsilon }_{n+1}\right)}^{+}||p_n|{|}^2.$$

$$\frac{1}{2}{\int }_{\Omega }|u_{n+1}{|}^2-|u_n{|}^2+|u_{n+1}-u_n{|}^{2}dx+{\int }_{\Omega }{k}_{n+1}\nu |\nabla u_{n+1}{|}^{2}dx$$

Full text

1 Introduction

Artificial compression (AC) methods are based on replacing $\nabla\cdot u=0$ by $\varepsilon p_{t}+\nabla\cdot u=0$ ($0<\varepsilon$ small), uncoupling velocity and pressure and advancing the pressure explicitly in time. Their high speed and low storage requirements recommend them for complexity bound fluid flow simulations. Unfortunately, time-accurate artificial compression approximations have proven elusive. Time accuracy (along with increased efficiency and decreased memory) is obtained by time-adaptive algorithms. To our knowledge, the defect correction based scheme of Guermond and Minev [17] and the non-autonomous AC method in [6], both adapting the time-step with $\varepsilon=k$ (time-step), are the only previous implicit, time-adaptive AC methods.

This report presents time-adaptive AC algorithms based on a new approach of independently adapting the AC parameter $\varepsilon$ and time-step $k$. The methods proceed as follows. A standard, first-order, implicit method, (1st Order) below, is used to advance the momentum equation in the artificial compression equations. A second-order velocity approximation, (2nd Order) below, is then computed at negligible cost using a time filter adapted from [19]. The difference between the first-order and second-order approximations gives a reliable estimator, EST(1), for the local error in the momentum equation for the first-order method and is used to adapt the time step in Algorithm 4.1, Section 4.

Adapting the AC parameter $\varepsilon$ is more challenging. Stability of the standard AC discrete continuity equation ($\varepsilon p_{t}+\nabla\cdot u=0$) is unknown for variable $\varepsilon$, [6]. We present two new, variable $\varepsilon,$ discrete continuity equations in (1.4) below and prove their unconditional, long-time stability in Theorems 2.1, 2.2 and 3.2. These results show that adaptivity will respond to accuracy constraints rather than try to correct stability problems with small time-steps. In these continuity equations, the size of $||\nabla\cdot u||$ is monitored and used to adapt the choice of the AC parameter $\varepsilon$ (e.g., Algorithm 3.1, Section 3) whereupon the calculation proceeds to the next time step. The self-adaptive strategy for independently adapting $\varepsilon$ also side steps the practical problem of how to pick $\varepsilon$ in AC methods and related penalty methods, even for constant time-steps. The new discrete continuity equations reduce to the standard $\varepsilon p_{t}+\nabla\cdot u=0$ for constant $\varepsilon$, improve, through greater simplicity, a non-autonomous ($\varepsilon=\varepsilon(t)$) AC formulation in [6] and yield now three proven stable extensions of the discrete AC continuity equation to variable $\varepsilon$. A comparison of the three is presented in Section 5. Determining if one or some combination of the three111The stability proof extends to weighted averages of the three discrete continuity equations. or some other, yet undetermined, possibility is to be preferred is an important open problem.

The second-order method. To obtain an $\mathcal{O}(k^{2})$ approximation of the momentum equation (with embedded error estimator), Algorithms 4.1 and 4.2 incorporate a recent idea of [19] of increasing accuracy and estimating errors by time filters. Theorem 3.2 of Section 3.1 gives a proof of unconditional, long-time stability of the second order, constant time-step but variable $\varepsilon$ method. The resulting embedded structure of Algorithms 4.1 and 4.2 suggests low-complexity, variable-order methods may be possible once an adaptive $\varepsilon$ strategy is well developed.

The second-order method is a one leg method. Reliable estimators of the local truncation error (LTE) in one leg methods are expensive as detailed in [10]. An inexpensive estimator, EST(2) in Algorithm 4.2, of the LTE in the method’s linear multistep twin, based on a second time filter, is presented. For the one leg method, this estimator is inexpensive but heuristic. The doubly adapted, second-order method in Algorithm 4.2 is tested in Section 5. The embedded structure of the first and second-order method suggests that adapting the method order in addition to the time-step and AC parameter $\varepsilon$ may increase accuracy and efficiency further.

Three stable treatments of the momentum equation (first, second and even variable order) are possible. Three stable treatments of the variable $\varepsilon$ continuity are now possible: two in (1.1) below and one in [6]. The result is nine adaptive AC methods with computational complexity comparable to the common first-order method, described next.

1.1 Review of a Common Artificial Compression Method

Denote by $u$ the velocity, $p$ the pressure, $\nu$ the kinematic viscosity, and $f$ the external force. Consider the slightly compressible/hyposonic222We do not include a traditional superscript ”$\varepsilon$” as we shall focus only on AC models and methods., [38], approximation to the incompressible Navier-Stokes equations in a domain $\Omega$ in ${{\mathbb{R}}}^{d},d=2,3$

[TABLE]

This is the most common of several possible formulations reviewed in Section 1.1 of [6]. To present methods herein we will consistently suppress the secondary spacial discretization333All stability results proven herein hold, by the same proof, for standard variational spatial discretizations such as finite element methods with div-stable elements.. Let $u^{\ast}$ denote the standard (second order) linear extrapolation of $u$ from previous values444Temperton and Staniforth [33] advocated even higher order extrapolation. to $t_{n+1}$

[TABLE]

To fix ideas, among many possible, e.g., [14], [15], [16], [22], [24], [27], [9], [26], [37], consider a common, constant time-step, semi-implicit time discretization of (1):

[TABLE]

Here $k$ is the time-step, $t_{n}=nk,$ $u_{n}$, $p_{n}$ are approximations to the velocity and pressure at $t=t_{n}$. This has consistency error $\mathcal{O}(k+\varepsilon)$ leading to the most common choice of selecting $\varepsilon=k$ to balance errors. Since $\nabla p_{n+1}=$ $\nabla p_{n}-(k/\varepsilon)\nabla\nabla\cdot u_{n+1}$, this uncouples into a velocity solve followed by an algebraic pressure update

[TABLE]

For constant $\varepsilon,k$, this method is unconditionally, nonlinearly, long-time stable, e.g., [14], [15], [31], [30]. Its long-time stability for variable $\varepsilon,k$ is an open problem, [6].

1.2 New Methods for Variable $\varepsilon,k$

Although well motivated, the  choice $\varepsilon=k$ cannot be more than a step to a correct choice. First observe that $Units(\varepsilon)=Time^{2}/Length^{3}$ while $Units(k)=Time$. Thus, a correct choice of $\varepsilon$ should be scaled to be dimensionally consistent and afterwards the constant multiplier optimized. Aside from dimensional inconsistency, the standard choice $\varepsilon=k$ ignores the different roles of $\varepsilon$ and $k$. To leading orders, the consistency error in the continuity equation is $\mathcal{O}(\varepsilon)$, independent of $k$, and the consistency error in the momentum equation is $\mathcal{O}(k)$, independent of $\varepsilon$. This observation on the standard method (2), (3) motivates the development plan for the doubly adaptive algorithms herein:

• •

Develop first (Section 2) and second (Section 3) order methods stable for variable $k,\varepsilon$.

• •

Adapt $\varepsilon_{n}$ to control the consistency error in the continuity equation by monitoring $||\nabla\cdot u||$, Sections 3, 4.

• •

Develop inexpensive estimators for momentum equation consistency error and adapt $k=k_{n}$ for its control, Section 4.

• •

Use (Section 4) and test (Section 5) the estimators in a doubly adaptive, variable $\varepsilon,k$, algorithm.

In adaptive methods, strong stability is necessary, so $\varepsilon_{n},k_{n}$ can be adapted for time-accuracy rather than to correct instabilities. One key difficulty, resolved by the two methods (4) below, is that useful stability is unknown for the common AC method (2) with variable $\varepsilon$, see [6], and even for the continuum model (1) with $\varepsilon=\varepsilon(t)$. A second key difficulty is that (unconditional, nonlinear) G-stability for variable time-steps is uncommon555To our knowledge, the only such two-step method is the little explored one of Dahlquist, Liniger, and Nevanlinna [7]. This second issue may be resolvable by a variable (first and second) order implementation since it would include the A-stable, fully implicit method.. (For example, the popular BDF2 method loses $A-$stability for increasing time-steps.)

The continuity equation is treated by either a geometric average (GA-Method) or a minimum term (min-Method) as follows. Given $u_{n},p_{n},\varepsilon_{n}$, select $\varepsilon_{n+1},k_{n+1}$ calculate $u_{n+1}$ then666A convex combination of the two continuity equations discretizations is also stable.

[TABLE]

These methods are proven in Section 2 to be unconditionally, variable $\varepsilon,k$ stable. For the discrete momentum equation, recall $u^{\ast}$ is an extrapolated approximation to $u(t_{n+1})$. The first-order method’s momentum equation is the standard one (2) above given by

[TABLE]

The (linearly implicit) treatment of the nonlinear term is inspired by Baker [4]. The second method, adapted from [19], adds a time filter to obtain $\mathcal{O}(k^{2})$ accuracy and automatic error estimation as follows. Let the time-step ratio be denoted $\tau=k_{n+1}/k_{n}$. Call $u_{n+1}^{1}$ the solution obtained from the first-order method (1st Order) above. The second-order approximation $u_{n+1}$ is obtained by filtering $u_{n+1}^{1}$:

[TABLE]

Denote by $D_{2}(n+1)$ the quantity above in braces

[TABLE]

Note that $D_{2}(n+1)$ is $2k_{n}k_{n+1}\times$(a second divided difference).

The usual $L^{2}$ norm $||\cdot||$ and inner product $(\cdot,\cdot)$ are denoted

[TABLE]

A simple estimate of the local error in the first-order approximation $u_{n+1}^{1}$ is given by a measure (here the $L^{2}$ norm) of the difference of the two approximations

[TABLE]

Estimating the error in the second-order approximation. Naturally one would like to use the second-order approximation for more than an estimator. It is possible to use $EST(1)$ above as a pessimistic estimator for $u_{n+1}$. In Section 3 we show that, eliminating the intermediate step $u_{n+1}^{1}$, the second-order method is equivalent to the second-order, one leg method (13) below. Estimation of the LTE for this OLM cannot be done by a simple time filter for reasons delineated in [10] and based on classical analysis of the LTE in OLMs of Dahlquist. We test an inexpensive but heuristic estimator that can be calculated by a second time filter. $EST(2)$ below is an LTE estimator for the OLMs linear multi-step twin. To estimate the local error in the second order approximation we use the third divided difference with multiplier chosen (by a lengthy but elementary Taylor series calculation) to cancel the first term of the LTE of the methods linear multi-step twin

[TABLE]

The resulting adaptive algorithm uncouples like (3) into a velocity update with a grad-div term then an algebraic pressure update. More reliable but more expensive estimators are possible. The above inexpensive but heuristic one is tested herein because the motivation for AC methods is often based on the need for faster and reduced memory algorithms in specific applications.

Section 2 presents the analysis of the two first-order methods, proving long-time, unconditional stability for variable $\varepsilon,k$. This analysis develops the key treatment of the discrete continuity equation necessary for stability. Section 3.1 gives a proof of unconditional, long time stability for the variable $\varepsilon$, constant $k$ second order method. This proof can be extended to decreasing time-steps but not increasing time-steps.

1.3 Related work

Artificial compression (AC) methods were introduced in the 1960’s by Chorin, Oskolkov and Temam. Their mathematical foundation has been extensively developed by Shen [29], [30], [31], [32] and Prohl [27]. Recent work includes [24], [9], [15], [16], [22], [26] and [37]. The GA-method (geometric averaging method) herein is motivated by work in [5] for uncoupling atmosphere-ocean problems stably.

There has been extensive development of adaptive methods for assured accuracy in fully coupled, $v-p$ discretizations, e.g., [21], and adaptive methods based on estimates of local truncation errors including [20], [23], [34]. In complement, the work herein aims at methods that use less expensive local (rather than global) error estimators, do not provide assured time-accuracy but emphasize (consistent with the artificial compression methods) low cognitive, computational, and space complexity. Aside from [6] and Guermond and Minev [17], extension of implicit, time-adaptive methods to artificial compression discretizations is undeveloped.

Herein accuracy is increased and local errors estimated by time filters. Other approaches are clearly possible. Time filters are an important tool in GFD to correct weak instabilities and extend forecast horizons, [3], [25], [28], [35], [36]. In [19], it was noticed that a time filter can also increase the convergence rate of the backward Euler method and estimate errors. G-stability of the resulting (constant time-step) time discretization was recently proven for the* fully-coupled*, velocity-pressure Navier-Stokes equations in [11].

2 First-Order, Variable $k,\varepsilon$ Methods

This section establishes unconditional, long-time, nonlinear stability of the two variable $k,\varepsilon$ first-order methods of Section 1.2 in the usual $L^{2}(\Omega)$ norm, denoted $||\cdot||$ with associated inner product $(\cdot,\cdot)$. The methods differ in the treatment of the discrete continuity equation and reduce to the standard AC method (2) for constant $\varepsilon,k$. We prove that the first order implicit discretization of the momentum equation with both new methods (5), (6) are unconditionally, nonlinearly, long-time stable without assumptions on $\varepsilon_{n},k_{n}$. We study these new methods in a bounded, regular domain $\Omega$ subject to the initial and boundary conditions

[TABLE]

The two, first-order methods are: Given $u_{n},p_{n},\varepsilon_{n}$, select $\varepsilon_{n+1},k_{n+1}$ and

[TABLE]

For constant $\varepsilon$ both methods reduce to the standard method (2), (3) for which stability is known. Thus, the interest is stability for variable $\varepsilon$.

Stability of the min-Method. It is useful to recall that

[TABLE]

Theorem 1** (Stability of the min Method).**

The variable $\varepsilon,k$ min-Method is unconditionally, long-time stable. For any $N>0$ the energy equality holds:

[TABLE]

Consequently, the stability bound holds:

[TABLE]

Proof.

First we note that using the polarization identity, algebraic rearrangement and considering the cases $\varepsilon_{n+1}>\varepsilon_{n}$ and $\varepsilon_{n+1}<\varepsilon_{n}$ we have

[TABLE]

We have $\varepsilon_{n+1}-\min\{\varepsilon_{n+1},\varepsilon_{n}\}=\left(\varepsilon_{n+1}-\varepsilon_{n}\right)^{+}$ and $\varepsilon_{n}-\min\{\varepsilon_{n+1},\varepsilon_{n}\}=$

$\left(\varepsilon_{n}-\varepsilon_{n+1}\right)^{+}.$ Thus,

[TABLE]

With this identity, take the inner product of the first equation with $k_{n+1}u_{n+1}$, the second with $k_{n+1}p_{n+1}$, integrate over the flow domain, integrate by parts, use skew symmetry, use the polarization identity twice and add. This yields

[TABLE]

From (7) the energy equality becomes

[TABLE]

Upon summation the first two terms telescope, completing the proof of the energy equality. The stability estimate follows from the energy equality and the Cauchy-Schwarz-Young inequality.

The stability analysis shows that the numerical dissipation in the min-Method is

[TABLE]

**The GA-Method. **The proof of stability of the GA-method differs from the last proof only in the treatment of the variable $\varepsilon$ term, resulting is a different numerical dissipation for the method.

Theorem 2** (Stability of GA-Method).**

The variable $\varepsilon,k$, first-order GA-Method is unconditionally, long-time stable. For any $N>0$ the energy equality holds:

[TABLE]

and the stability bound holds:

[TABLE]

Proof.

First we note that using the polarization identity we have

[TABLE]

The remainder of the proof is the same as for the min-Method.

The stability analysis shows that the numerical dissipation in the GA-Method is

[TABLE]

There is no obvious way to tell á priori which method’s numerical dissipation is larger or to be preferred. A numerical comparison is thus presented in Section 5.

Remark 3**.**

The continuum analogs. It is natural to ask if there is a non-autonomous continuum AC model associated with each method. The momentum equation for each continuum model is the standard

[TABLE]

The associated continuum continuity equation for the min-Method is

[TABLE]

whereas the continuum continuity equation for the GA-method is

[TABLE]

Analyzing convergence of each to a weak solution of the incompressible NSE as (non-autonomous) $\varepsilon(t)\rightarrow 0$ is a significant open problem.

3 Second-Order, Variable $\varepsilon$ Methods

The first-order methods are now extended to embedded first and second-order methods adapting [19] from ODEs to the NSE. First we review the idea of extension used.

Review of the ODE algorithm. Consider the initial value problem

[TABLE]

Recall $\tau=k_{n+1}/k_{n}$ is the time-step ratio. The second-order accurate, variable time-step method of [19] is the standard backward Euler (fully implicit) method followed by a time filter:

[TABLE]

The combination is second-order accurate, $A-$stable for constant or decreasing time-steps and a measure of the pre- and post-filter difference

[TABLE]

can be used in a standard way as a local error estimator for the lower order approximation $y_{n+1}^{1}$ or a (pessimistic) estimator for the higher order approximation $y_{n+1}$.

A simple, adaptive$-\varepsilon$, second-order AC algorithm. The continuity equation for both methods can be written

[TABLE]

This can be used to uncouple velocity and pressure using

[TABLE]

The discrete momentum equation for either first-order method is then

[TABLE]

Applying the time filter of (10) to the velocity approximation increases the methods accuracy to $\mathcal{O}(k^{2})$. This combination yields a simple, second-order, constant time-step but adaptive $\varepsilon$ algorithm. In the algorithm below the change in $\varepsilon$ is restricted to be between halving and doubling the previous $\varepsilon$ value.

Algorithm 4**.**

**[Simple, adaptive $\varepsilon$, constant time-step, second-order AC method]. Given $u_{n},u_{n-1},p_{n},k,\varepsilon_{n+1},\varepsilon_{n}$, and tolerance $TOL_{c}$, **

*Select:  $\hat{\varepsilon}=\sqrt{\varepsilon_{n+1}\varepsilon_{n}}$ $\ \mathtt{or}$ $\hat{\varepsilon}=\min\{\varepsilon_{n+1},\varepsilon_{n}\}$  *

*Set:  $u^{\ast}=2u_{n}-u_{n-1}.$  *

Solve for $u_{n+1}^{1}$

[TABLE]

Filter, Compute estimator $EST_{c}$ , Find $p_{n+1}$

[TABLE]

Adapt $\varepsilon:$  IF $EST_{c}>TOL_{c}$ , THEN repeat step after resetting $\varepsilon_{n+1}$ by

[TABLE]

ELSE

[TABLE]

and proceed to next step.

3.1 Stability of the second-order method for variable $\varepsilon$, constant $k$

This section establishes unconditional, nonlinear, long-time stability of the second-order GA-method for constant time-steps but variable $\varepsilon$. The proof addresses the interaction between the filter step with the continuity equation. It is adapted to the min-Method following ideas in the proof of Theorem 2.1. For constant time-steps and variable $\varepsilon$ the GA-method is as follows. Given $u_{n},p_{n},\varepsilon_{n}$, select $\varepsilon_{n+1}$ and $u^{\ast}=$ $2u_{n}-u_{n-1}$ (since the time-step is here constant). Then,

[TABLE]

We now prove an energy equality for the method which implies stability.

Theorem 5**.**

The method (12) satisfies the following discrete energy equality (from which stability follows). For any $N>1$

[TABLE]

Proof.

To prove stability, eliminate the intermediate value $u_{n+1}^{1}$ in the momentum equation. From the filter step $u_{n+1}=u_{n+1}^{1}-\frac{1}{3}\left\{u_{n+1}^{1}-2u_{n}+u_{n-1}\right\}$ we have

[TABLE]

Replacing $u_{n+1}^{1}$ by $\frac{3}{2}u_{n+1}-u_{n}+\frac{1}{2}u_{n-1}$ yields the equivalent discrete momentum equation:

[TABLE]

Multiply by the time-step $k$, take the $L^{2}$ inner product of the momentum equation (13) with $\frac{3}{2}u_{n+1}-u_{n}+\frac{1}{2}u_{n-1}$, the $L^{2}$ inner product of the discrete continuity equation with $p_{n+1}$ and add. Two pressure terms cancel since $u_{n+1}^{1}=\frac{3}{2}u_{n+1}-u_{n}+\frac{1}{2}u_{n-1}$and the nonlinear terms vanish due to skew-symmetry. Thus, we obtain

[TABLE]

The key terms are the first two. For the first term, apply the following identity from [11] with $a=u_{n+1},$ $b=u_{n},$ $c=u_{n-1}$

[TABLE]

This yields

[TABLE]

For the pressure term $\left(\sqrt{\varepsilon_{n+1}\varepsilon_{n}}p_{n},p_{n+1}\right)$ the polarization identity, suitably applied, yields

[TABLE]

Thus

[TABLE]

Combining the pressure and velocity identities, we have

[TABLE]

Summing from $n=1$ to $N$ proves unconditional, long-time stability.

4 Doubly $k,\varepsilon$ Adaptive Algorithms

We present three doubly adaptive AC algorithms: first-order, second-order method and a third that adapts the method order. The first two are tested in Section 5. While not tested herein, we include the variable order adaptive algorithm for its clear interest. In the first algorithm, the error is estimated by a time filter and the next time-step and next $\varepsilon$ are adapted777The formula for $\varepsilon_{new}$ could be improvable. based on

[TABLE]

In our implementation, a safety factor of $0.9$ is used and the maximum change in both is (additionally) restricted to be between $0.5$ & $2.0$.

Algorithm 6** **(Doubly $k$, $\varepsilon$ Adaptive, First-Order

Method).

*Given $TOL_{m}$,

$TOL_{c}$, $u_{n}$, $u_{n-1}$, $u_{n-2}$ and $k_{n+1}$, $k_{n}$, $k_{n-1}$ *

Compute: $\tau=\frac{k_{n+1}}{k_{n}}$ and  $\alpha_{1}=\frac{{\LARGE\tau(1.0+\tau)}}{{\LARGE 1.0+2.0\tau}}$

*Select:  $\hat{\varepsilon}=\sqrt{\varepsilon_{n+1}\varepsilon_{n}}$ $\mathtt{or}$ $\hat{\varepsilon}=\min\{\varepsilon_{n+1},\varepsilon_{n}\}.$  *

Set   $u^{\ast}=\left(1+\tau\right)u_{n}-\tau u_{n-1}.$

Find BE approximation $u_{n+1}$

[TABLE]

Compute difference $D_{2}$ and Estimators

[TABLE]

IF $EST_{c}>TOL_{c}$  or $EST(1)>TOL_{m}$ THEN repeat step after resetting $\varepsilon_{n+1},k_{n+1}$ by

[TABLE]

ELSE Predict best next step for each approximation:

[TABLE]

ENDIF

Update pressure:   $p_{n+1}=\frac{\hat{\varepsilon}}{\varepsilon_{n+1}}p_{n}-\frac{k_{n+1}}{\varepsilon_{n+1}}\nabla\cdot u_{n+1}.$

Proceed to next step.

The second-order, doubly adaptive algorithm. For the second-order, doubly adaptive method, we predict the next $\varepsilon$ value the same as in the first-order method and predict the next time step based on

[TABLE]

$EST(2)$ is calculated as follows. The second-order method is equivalent, after elimination of the intermediate (first-order) approximation, to a one leg method exactly as in (12) in the constant time-step case. The one leg method’s linear multistep twin has local error proportionate to $k^{3}u_{ttt}+O(k^{4})$. Thus, an estimate of $u_{ttt}$ is computed using difference of $D_{2}$ as follows. Write

[TABLE]

From differences of $D_{2}(n+1)$, $D_{2}(n)$ we obtain the estimator:

[TABLE]

where the coefficient $\alpha_{2}$ is determined through a Taylor series calculation to be

[TABLE]

Algorithm 7** (Doubly Adaptive, Second-Order Algorithm).**

*Given $TOL_{m}$, $TOL_{c}$, $u_{n}$, $u_{n-1}$, $u_{n-2}$, previous 2nd difference $D_{2}(n)$ and $k_{n+1}$, $k_{n}$, $k_{n-1}$ *

Compute: $\tau=\frac{k_{n+1}}{k_{n}},$ $\alpha_{1}=\frac{{\LARGE\tau(1.0+\tau)}}{{\LARGE 1.0+2.0\tau}}$, $\alpha_{2}=\frac{\tau_{n}(\tau_{n+1}\tau_{n}+\tau_{n}+1)(4\tau_{n+1}^{3}+5\tau_{n+1}^{2}+\tau_{n+1})}{3(\tau_{n}\tau_{n+1}^{2}+4\tau_{n}\tau_{n+1}+2\tau_{n+1}+\tau_{n}+1)}$

Select:   $\hat{\varepsilon}=\sqrt{\varepsilon_{n+1}\varepsilon_{n}}$ $\mathtt{or}$ $\hat{\varepsilon}=\min\{\varepsilon_{n+1},\varepsilon_{n}\}.$

Set:   $u^{\ast}=\left(1+\tau\right)u_{n}-\tau u_{n-1}.$

Find BE approximation $u_{n+1}^{1}$

[TABLE]

Compute difference $D_{2}$ and update velocity

[TABLE]

Compute estimators

[TABLE]

IF $EST_{c}>TOL_{c}$  or $EST(2)>TOL_{m}$ THEN repeat step after resetting $\varepsilon_{n+1},k_{n+1}$ by

[TABLE]

ELSE Predict best next step:

[TABLE]

Update pressure:   $p_{n+1}=\frac{\hat{\varepsilon}}{\varepsilon_{n+1}}p_{n}-\frac{k_{n+1}}{\varepsilon_{n+1}}\nabla\cdot u_{n+1}.$

Proceed to next step.

The adaptive order, time-step and $\varepsilon$ algorithm. To adapt $\varepsilon,k$ and the method order we use the local truncation error indicators for the momentum and continuity equations, respectively,

[TABLE]

The algorithm computes two velocity approximations. The first $u^{1}$ is first-order and $A-$stable for all combinations of time-step and $\varepsilon$. The second $u$ is second-order, A-stable for constant (or decreasing) time-step but only $0-$stable for increasing time-steps. Variable (1 or 2) order is introduced as follows. The local error in each approximation is estimated. If both are above the tolerance, the step is repeated. Otherwise, the optimal next time-step is predicted for each method by

[TABLE]

The actual $k_{n+1}$ presented below and in the tests in Section 5 is restricted to be $(0.5$ to $2.0)\times k_{n}$ and includes a safety factor of $0.9$.

Algorithm 8** (Adaptive order, $k$, $\varepsilon$).**

*Given $TOL_{m}$, $TOL_{c}$, $u_{n}\mathtt{,\ }u_{n-1}$, $u_{n-2}$, previous second difference $D_{2}(n)$ and $k_{n+1}$, $k_{n}$, $k_{n-1}$ *

Compute: $\tau=\frac{k_{n+1}}{k_{n}},$ $\alpha_{1}=\frac{{\LARGE\tau(1.0+\tau)}}{{\LARGE 1.0+2.0\tau}},$ $\alpha_{2}=\frac{\tau_{n}(\tau_{n+1}\tau_{n}+\tau_{n}+1)(4\tau_{n+1}^{3}+5\tau_{n+1}^{2}+\tau_{n+1})}{3(\tau_{n}\tau_{n+1}^{2}+4\tau_{n}\tau_{n+1}+2\tau_{n+1}+\tau_{n}+1)}$

Select:   $\hat{\varepsilon}=\sqrt{\varepsilon_{n+1}\varepsilon_{n}}$ $\mathtt{or}$ $\hat{\varepsilon}=\min\{\varepsilon_{n+1},\varepsilon_{n}\}.$

Set:   $u^{\ast}=\left(1+\tau\right)u_{n}-\tau u_{n-1}.$

Find BE approximation $u_{n+1}^{1}$

[TABLE]

Compute difference $D_{2}$ and updated velocity

[TABLE]

Compute estimators

[TABLE]

IF $EST_{c}>TOL_{c}$  or $\min\{EST(1),EST(2)\}>TOL_{m}$ THEN repeat step, resetting $\varepsilon_{n+1},k_{n+1}$ by

[TABLE]

ELSE Predict $\varepsilon,k$ for each approximation:

[TABLE]

Select method order with larger next step:

IF ($STEPBE>STEPFilter$) Then

$\ \ \ \ \ \ \ \ u_{n+1}=$*  $u_{n+1}^{1}$*

$k_{n+2}=STEPBE$**

ELSE      $\ k_{n+2}=STEPFilter$

ENDIF

Update pressure:  $p_{n+1}=\frac{\hat{\varepsilon}}{\varepsilon_{n+1}}p_{n}-\frac{k_{n+1}}{\varepsilon_{n+1}}\nabla\cdot u_{n+1}.$

Proceed to next step

The fixed order methods can, if desired, be implemented by commenting out parts of the variable order Algorithm 4.3.

5 Three Numerical Tests

The stability and accuracy of the new methods are interrogated in two numerical tests and the three discrete continuity equations are compared in our third test. The tests employ the finite element method to discretize space, with Taylor-Hood ($\mathbb{P}_{2}/\mathbb{P}_{1}$) elements, [18]. All the stability results proven herein hold for this spatial discretization by essentially the same proofs. The meshes used for both tests are generated using a Delaunay triangulation. The software package FEniCS is used for both experiments [1].

We begin with comparative tests of the adaptive $k,\varepsilon$, first and second-order method. Both adapt $\varepsilon$ based on $||\nabla\cdot u||$. The first-order method accepts the first-order approximation $u_{n+1}^{1}$ and adapts the time-step based on $EST(1)$. The second-order method accepts $u_{n+1}$ as the approximation and adapts the time step based on $EST(2)$.

5.1 Test 1: Flow Between Offset Circles

To interrogate stability and accuracy of the GA-method, we present the results of two numerical tests. Pick

[TABLE]

with no-slip boundary conditions on both circles and $\nu=0.001$. The finite element discretization has a maximal mesh width of $h_{max}=0.0133,$ and the flow was solved using the direct solver UMFPACK [8]. For this test, we use fixed tolerances $TOL_{m}=TOL_{c}=0.001$. The flow (inspired by the extensive work on variants of Couette flow, [12]), driven by a counterclockwise force (with $f\equiv 0$ at the outer circle), rotates about $(0,0)$ and interacts with the immersed circle. This induces a von Kármán vortex street which re-interacts with the immersed circle creating more complex structures. There is also a central (polar) vortex that alternately self-organizes then breaks down. Each of these events includes a significant pressure response.

For both approximations we track the evolution of $k_{n}$ and $\varepsilon_{n}$, the pressure at the origin, the violation of incompressibility, and the algorithmic energy $\|u_{h}^{n+1}\|^{2}+\varepsilon_{n+1}\|p_{h}^{n+1}\|^{2}$. These are all depicted in Figure 1 below.

Figure 1A shows that the second-order scheme consistently chooses larger time-steps than the first-order method. The evolution of $\varepsilon$, in Figure 1B, behaves similarly for both methods once the flow evolves. In testing AC methods pressure initialization often causes irregular, transient spiky behavior near $t=0$ such as in Figures 1A, 1B, 1D.

The behavior of the pressure at the origin, $p(0,0;t)$ vs. $t$, is depicted in Figure 1C. To our knowledge, there is no convergence theory for AC methods (or even fully coupled methods) which implies maximum norm convergence for the pressure over significant time intervals and for larger Reynolds numbers. Still, the irregular behavior observed in approximate solutions, while not conforming to a convergence theory, reflects vortex events across the whole domain and is interesting to compare. The profiles of the pressure at the origin are similar for both methods over $0\leq t\leq 4$. For $t>4,$ $p(0,0;t)$ for the second-order scheme is less oscillatory. This is surprising because the first-order scheme has more numerical dissipation. The divergence evolution of the schemes also differ in the initial transient of $||\nabla\cdot u(t)||$. After the initial transient, the divergence behavior is similar. It is also possible that the difference in $||\nabla\cdot u||$ transients is due to the strategy of $\varepsilon-$adaptation being sub-optimal. The model energy of both methods is largely comparable. We note that the model energy depends on the choices of $\varepsilon$ made. Thus model energy is not expected to coincide exactly. Generally, Figures 1D–1E behave similarly for both algorithms.

5.2 Test 2: Convergence and Adaptivity

The second numerical test concerns the accuracy and adaptivity of the GA-method. Let $\Omega=]0,1[^{2}$, with $\nu=1$. Consider the exact solution (obtained from [14] and applied to the Navier-Stokes equations)

[TABLE]

and consider a discretization of $\Omega$ obtained by 300 nodes on each edge of the square. We proceed by running five experiments, adapting both the first- and second-order schemes using the algorithms above, where the tolerance for the continuity and momentum equations is $10^{-(.25i+3)}$ for $i=0,1,2,3,4$. To control the size of the timesteps, we require $k_{n}$ to be chosen such that $EST(1)\in(TOL_{m}/10,TOL_{m})$. The solutions were obtained in parallel, utilizing the MUMPS direct solver [2]. To examine convergence, we present in Figure 2 log-log plots of the errors of the pressure and the velocity against the average time-step taken during the test. We also present semilog plots of the evolution of the pressure error and timestep during the final test below. The plots show that the time-step adaptation is working as expected and reducing the velocity error, Figure 2C. Our intuition is that the pressure error is linked to satisfaction of incompressibility; however, Figure 2D indicates convergence with respect to the timestep. In our calculations we did observe the following: If $||\nabla\cdot u||$ is, e.g., two orders of magnitude smaller then the tolerance, $\varepsilon$ is rapidly increased to be even $\mathcal{O}(1)$. At this point the pressure error and violation of incompressibility spike upward and $\varepsilon$ is then cut rapidly. This behavior suggests that a band of acceptable $\varepsilon$-values should be imposed in the adaptive algorithm.

To compare the GA, Min method and the scheme introduced in [6], we use the test problem given above in this section with a known exact solution. The results are given in Figure 3 below. Here, we use a mesh with the same density and final time $T=1$. A timestep $k_{n}=10^{-2}$ is kept constant in this run to highlight differences in the evolution of the variable $\varepsilon_{n}$, which has an initial value $\varepsilon_{0}=10^{-4}$. These tests are preliminary: In them, the min-Method seems preferable in error behavior but yields smaller values and thus less well-conditioned systems. In the evolution of all four quantities, the GA- and the CLM [6] method exhibit near identical behavior. The min-Method, however, forces $\varepsilon$ to be an order of magnitude lower than the values obtained by the other two schemes. This, in turn, forces the divergence to be reduced. Furthermore, both the velocity and pressure errors for the min-Method are smaller than those of the GA- and CLM-Methods.

6 Conclusions, open problems and future prospects

There are many open problems and algorithmic improvements possible. The doubly adaptive algorithm selected smaller values of $\varepsilon$ than $k$ in our tests with the same tolerance for both. A further synthesis of the methods herein with the modular grad-div algorithm of [13] would eliminate any conditioning issues in the linear system arising. Developing doubly adaptive methods of order greater than two (with modular grad-div) is an important step to greater time accuracy. We mention in particular the new embedded family of orders 2,3,4 of [10] as a natural extension. The method of Dahlquist, Liniger and Nevanlinna [7] is unexplored for PDEs, but has promise in CFD because it is A-stable for both increasing and decreasing time-steps. Improved error estimators for the second-order method herein would increase reliability. For AC methods, pressure initialization and damping of nonphysical acoustics are important problems where further progress would be useful.

Open problems. The idea of adapting independently $k$ and $\varepsilon$ is promising but new so there are many open problems. These include:

• •

Is the $\varepsilon$-adaptation formula $\varepsilon_{new}=\varepsilon_{old}(TOL/||\nabla\cdot u||)$ improvable? Perhaps the quotient should be to some fractional power. Perhaps adapting $\varepsilon$ should be based of a relative error in $||\nabla\cdot u||$, such as $||\nabla\cdot u||/||\nabla u||$. Analysis of the local (in time) error in $||\nabla\cdot u||$ is needed to support an improvement.

• •

The $\varepsilon$-adaptation strategy seems to need preset limits, $\varepsilon_{\min},\varepsilon_{\max}$, to enforce $\varepsilon_{\min}\leq\varepsilon\leq\varepsilon_{\max}$. The preset of $\varepsilon_{\min}$ is needed because $\nabla\cdot u=0$ cannot be enforced pointwise in many finite element spaces. Finding a reasonable strategy for these presets is an open problem. Similarly, it would be useful to develop a coherent strategy for relating the two tolerances rather than simply picking them to be equal (as herein).

• •

Proving convergence to a weak solution of the incompressible NSE of solutions to the continuum analogs of the GA-method and min-Method for variable $\varepsilon$ is an important open problem. In this analysis it is generally assumed that $\varepsilon(t)\rightarrow 0$ in an arbitrary fashion. A more interesting problem is to link $\varepsilon(t)$ and $||\nabla\cdot u||$ in the analysis. Similarly, an á priori error analysis for variable $\varepsilon$ is an open problem and may yield insights on how the variance of $\varepsilon(t)$ should be controlled within an adaptive algorithm. The consistency error of the two methods are $\mathcal{O}(k+\varepsilon)$ and $\mathcal{O}(k^{2}+\varepsilon)$, respectively. Energy stability has been proven herein for the first order method and for the constant time-step, second order method. Thus, error estimation while technical, should be achievable.

• •

Comprehensive testing of the variable (first or second) order method is an open problem. VSVO methods are the most effective for systems of ODEs but have little penetration in CFD. Testing the relative costs and accuracy of VSVO in CFD is an important problem.