# Improved convergence of the Arrow-Hurwicz iteration for the Navier-Stokes
equation via grad-div stabilization and Anderson acceleration

Pelin G. Geredeli<EMAIL_ADDRESS>Leo G. Rebholz<EMAIL_ADDRESS>Duygu Vargun<EMAIL_ADDRESS>Department of Mathematical Sciences, Clemson
University, Clemson, SC 29634, USA  Department of Mathematics, Iowa State
University, Ames IA, 50011, USA  Ahmed Zytoon<EMAIL_ADDRESS>
###### Abstract

We consider two modifications of the Arrow-Hurwicz (AH) iteration for solving
the incompressible steady Navier-Stokes equations for the purpose of
accelerating the algorithm: grad-div stabilization, and Anderson acceleration.
AH is a classical iteration for general saddle point linear systems and it was
later extended to Navier-Stokes iterations in the 1970’s which has recently
come under study again. We apply recently developed ideas for grad-div
stabilization and divergence-free finite element methods along with Anderson
acceleration of fixed point iterations to AH in order to improve its
convergence. Analytical and numerical results show that each of these methods
improves AH convergence, but the combination of them yields an efficient and
effective method that is competitive with more commonly used solvers.

###### keywords:

Anderson acceleration, Arrow-Hurwicz, Navier-Stokes equations, Finite element
method (FEM)

##  1 Introduction

We consider in this paper solving the incompressible steady Navier-Stokes
equations (NSE) with the Arrow-Hurwicz (AH) iteration. The steady NSE defined
on a domain $\Omega\subset\mathbb{R}^{d}$ ($d$=2 or 3) are given by

$\displaystyle-\nu\Delta{u}+(u\cdot\nabla)u+{\nabla}p$ $\displaystyle=f\qquad$
$\displaystyle\text{in }\Omega,$ (1.1a) $\displaystyle\nabla\cdot u$
$\displaystyle=0\qquad$ $\displaystyle\text{in }\Omega,$ (1.1b) $\displaystyle
u$ $\displaystyle=0\qquad$ $\displaystyle\text{on }{\partial}\Omega,$ (1.1c)

where $u$ and $p$ represent the unknown velocity and pressure, $f$ a given
forcing, and $\nu$ the kinematic viscosity which is inversely proportional to
the Reynolds number $Re$.  Extension of this work to one time step in a
temporal discretization of the time dependent NSE is straight-forward.

Among various novel iterative methods for solving saddle point systems, the AH
algorithm for the steady NSE was seemingly first studied by Temam in 1977 in
[37], and also  more recently in [7]. The AH iteration is given with the
following decoupled equations:

$\displaystyle-\frac{1}{\rho}\Delta(u^{m+1}-u^{m})-\nu\Delta
u^{m}+u^{m}\cdot\nabla u^{m+1}+\nabla p^{m}$ $\displaystyle=f,$
$\displaystyle\alpha(p^{m+1}-p^{m})+\rho\nabla\cdot u^{m+1}$
$\displaystyle=0,$

with $\rho$ and $\alpha$ being user determined parameters. We note that if
$\rho=\nu^{-1}$ then AH is exactly the modified Uzawa algorithm from [8]. This
iteration is interesting because it is efficient since the two equations
decouple, with the second equation being simple and the first equation
requiring a typical convection diffusion solver where one controls the
diffusion coefficient (in each iteration) with $\rho$. Hence from an
implementation (i.e. linear algebraic) point of view, the cost of one AH
iteration is very cheap compared to that of a typical Picard or Newton
iteration which needs to resolve a saddle point system. However, a serious
drawback of the AH method is that its convergence properties are not
particularly good and even though each iteration is cheap, the total number of
iterations can be very large. The purpose of this paper is to improve the AH
algorithm so that it is a competitive and even attractive method for efficient
computing of accurate steady NSE solutions. We enhance the AH method with two
recently developed ideas, one from computational fluid dynamics (grad-div
stabilization) and the other from nonlinear solver theory (Anderson
acceleration). Indeed, we show that the combination of these two improvements
theoretically and computationally yields that AH method can become a very good
solver.

The classical AH iteration developed in 1958 by Arrow and Hurwicz [3] is a
stationary iterative method to solve saddle point linear systems. As noted in
[4], this linear algebraic AH iteration can be regarded as an inexpensive
alternative to the (linear algebraic) Uzawa method [39] whenever solves with
the matrix arising from the convection-diffusion operators are expensive.
Temam seems to be the first to export the AH iteration ideas to Galerkin
methods for solving the steady NSE, and was able to prove convergence
(although without a rate) under certain choices of parameters. The iteration
was (finally) proven to be contractive in 2017 [7], where it was shown to be
linearly convergent under very small data and certain choices of parameters.
While this was a big step forward for AH and the convergent rate was proven
less than 1, the exact rate was not easy to decipher and in practice could be
very close to 1\. The numerical tests in [7] for some relatively easy problem
revealed that the AH method could be made to converge with good parameter
choices, but the number of iterations could be very large (e.g. over 700
iterations with $Re=100$ for a 2d driven cavity problem, and over 10,000 for a
2d steady flow past a cylinder). While iterations of AH would likely be 5-20
times cheaper than one iteration of usual Picard (e.g. if Krylov solvers with
preconditioners such as those in [19, 6, 5] to solve the saddle point linear
systems at each iteration), such high iteration counts still make AH
uncompetitive.

Herein we aim to improve the convergence properties of the AH method by
enhancing it with two techniques. The first is the addition of grad-div
stabilization  which refers to consistent penalization term that adds
$0=-\gamma\nabla(\nabla\cdot u)$ to the NSE momentum equation before
discretizing, where $\gamma>0$ is a user defined parameter (how large it
should be depends on many factors, see e.g. [22]). It was first proposed by
Hughes and Franca in 1988 [11], and has been shown to improve accuracy of
finite element approximations [28], improve saddle point linear solvers [28,
5, 19], and help various NSE nonlinear iterative solvers converge faster, e.g.
[35, 34]. The grad-div stabilized AH iteration takes the form

$\displaystyle-\frac{1}{\rho}\Delta(u^{m+1}-u^{m})-\nu\Delta
u^{m}+u^{m}\cdot\nabla u^{m+1}-\gamma\nabla(\nabla\cdot u^{m+1})+\nabla p^{m}$
$\displaystyle=f,$ $\displaystyle\alpha(p^{m+1}-p^{m})+\rho\nabla\cdot
u^{m+1}$ $\displaystyle=0.$

We show that the existing convergence theory can dramatically be improved with
the use of grad-div stabilization theory. Moreover, we show that under a
certain choice of parameters and in a particular (but commonly used) discrete
setting, grad-div stabilized AH method is equivalent to the classical iterated
penalty Picard iteration. By establishing this connection, we are able to
bring to bear the long established theory for this classical iteration to the
AH setting, which establishes a linear convergence rate close to that of
Picard for sufficiently large $\gamma$.

The second enhancement we provide to the AH method is Anderson acceleration
(AA). AA is an extrapolation technique used to improve convergence of fixed
point iterations. It was first developed in 1965 by D.G. Anderson [2], and its
use has exploded in the last decade after the paper of Walker and Ni in 2011
showed how effective AA can be on a wide range of problems [40]. It has
recently been used to improve convergence and robustness of solvers for
various types of flow problems [26, 31, 32], geometry optimization [29],
radiation diffusion and nuclear physics [1, 38], molecular interaction [36],
and many others e.g. [40, 23, 25, 26, 12, 41, 20]. In [33], AA was shown to
significantly improve the convergence and robustness for the IPP method for
the NSE and allow for a much wider range of penalty parameter choices. Given
the success AA has had in improving other types of nonlinear iterations for
the NSE, applying it to the AH method seems a natural next step. Moreover, due
to its dramatic improvement of the IPP method in [33] and our showing the
strong connection of grad-div stabilized AH method to IPP method, applying AA
to grad-div stabilized AH seems an optimal combination to improve AH
convergence behavior. A general convergence framework was developed for AA in
[10] and then sharpened in [30] which allows for theoretical justification of
improved linear convergence from AA, if the associated fixed point function
satisfies sufficient smoothness properties. We will set up the AH iteration as
a fixed point problem and prove that its fixed point operator satisfies the
assumptions needed to apply the AA convergence theory. Furthermore, extensive
computations of AH with AA are performed, and AA is observed to provide a
dramatic improvement in convergence behavior, with the best convergence coming
from combining AA with grad-div stabilization.

This article is arranged as follows. Section 2 provides the necessary notation
and mathematical preliminaries used throughout the paper. In section 3 we
consider the theoretical improvement provided by grad-div stabilization, while
in section 4 we show how the fixed point operator associated with the AH
iteration allows for the AA theory from [30] to be applied. Finally, in
section 5, we give results of several numerical tests that show AH method
enhanced with AA and grad-div stabilization can be a very effective nonlinear
solver for the steady NSE.

##  2 Preliminaries

In this section we provide some mathematical preliminaries and notation that
will be used throughout the paper. We begin by defining the following function
spaces on a domain $\Omega$ that either has smooth boundary or is a convex
polygon:

$\displaystyle L^{2}(\Omega)$ $\displaystyle:=\\{w:\Omega\mapsto\mathbb{R}:\
\|w\|_{{\it L}^{2}(\Omega)}:=\left(\int_{\Omega}\
{|w|^{2}}\,dx\right)^{1/2}<\infty\\},$ $\displaystyle H^{m}(\Omega)$
$\displaystyle:=\\{w:\Omega\mapsto\mathbb{R}:\ \|w\|_{{\it
H}^{m}(\Omega)}:=\left(\sum_{|\beta|\leq m}\|{\it D}^{\beta}w\|_{{\it
L}^{2}(\Omega)}^{2}\right)^{1/2}<\infty\\}.$

Throughout this paper, $(.,.)$ and $\|\cdot\|$ denote the inner product and
norm on $L^{2}({\Omega})$, respectively. All other norms will be denoted with
subscripts. Also, we define the following natural spaces for NSE:

$\displaystyle Q$ $\displaystyle:=\\{w\in{\it L}^{2}(\Omega):\ \int_{\Omega}\
{w}\,dx=0\\},$ $\displaystyle X$ $\displaystyle:=\\{{w}\in{\it
H}^{1}(\Omega):\ {w}|_{\partial\Omega}=0\\},$

and we do not distinguish vector and scalar valued spaces, as it will be clear
from context.

The skew-symmetric trilinear form $b^{*}:X\times X\times
X\rightarrow\mathbb{R}$ is defined by

$b^{*}(u,v,w)=\frac{1}{2}(((u\cdot\nabla)v,w)-(u\cdot\nabla)w,v)),$ (2.1)

and it can easily be observed that

$b^{*}(u,v,v)=0.$ (2.2)

We will also utilize the well known bound resulting from Hölder’s and Sobolev
inequalties [24]:

$\left|b^{*}(u,v,w)\right|\leq M\|\nabla u\|\|\nabla v\|\|\nabla w\|,$ (2.3)

where $M$ is a constant depending only on $\Omega$.

###  2.1 Finite element preliminaries

Let $X_{h}\times Q_{h}\subset X\times Q$ be conforming and finite dimensional
finite element spaces for the velocity and pressure. Then a finite element
method for (1.1), based on the standard velocity-pressure formulation and
equipped with grad-div stabilization seeks $(u,p)\in X_{h}\times Q_{h}$ such
that $\forall(v,q)\in X_{h}\times Q_{h}$ we have

$\displaystyle\nu(\nabla{u},\nabla{v})+\gamma(\nabla\cdot u,\nabla\cdot
v)+b^{*}(u,u,v)-(\nabla\cdot v,p)$ $\displaystyle=({f},v)$ (2.4a)
$\displaystyle(\nabla\cdot u,q)$ $\displaystyle=0.$ (2.4b)

Here, $\gamma\geq 0$ is referred to as the grad-div stabilization parameter.
The discrete problem (2.4) is well-posed if the pair satisfies the inf-sup
condition

$\sup_{0\neq v\in X_{h}}\frac{(\nabla\cdot v,q)}{\|\nabla
v\|}\geq\beta\|q\|\qquad\forall q\in Q_{h},$

for some $\beta>0$, and the small data condition
$\kappa:=M\nu^{-2}\|f\|{-1}<1$ holds. The small data condition is needed for
uniqueness (although precisely how sharp it is remains an open question),
although existence and boundedness can be proven for any given data. Common
choices that satisfy the inf-sup condition with $\beta$ independent of $h$,
and the ones we make herein are $X_{h}\times Q_{h}=P_{k}(\tau_{h})\cap X\times
P_{k-1}(\tau_{h})\cap C^{0}(\Omega)$ Taylor Hood elements (with $\tau_{h}$
representing a regular conforming mesh of $\Omega$) and $X_{h}\times
Q_{h}=P_{k}(\tau_{h})\cap X\times P_{k-1}(\tau_{h})\cap Q$ Scott-Vogelius
elements with appropriate $k$ and mesh structure (see e.g. [18, 22] for more
details).

We will utilize the bound

$\|\nabla u\|\leq\nu^{-1}\|f\|_{-1},$ (2.5)

which is proven in [14, 24] and can be easily deduced from (2.4).

###  2.2 The Arrow-Hurwicz method

We recall the Arrow-Hurwicz (AH) method from [7] for steady Navier-Stokes
equations. The method is given in [37] with a slight change of parameter
variables.

######  Algorithm 2.1.

Let $\rho,\alpha>0$ be user selected parameters, then

  1. 1.

Let $u^{0}\in X_{h}$ and $p^{0}\in Q_{h}$ be the solution of the mixed
formulation: $\forall({v},q)\in X_{h}\times Q_{h}$ :

$\displaystyle(\nabla{u^{0}},\nabla{v})-(\nabla\cdot v,p^{0})$
$\displaystyle=({f},v)$ (2.6a) $\displaystyle(\nabla\cdot u^{0},q)$
$\displaystyle=0,$ (2.6b)

  2. 2.

For $m\geq 0$, we define $u^{m+1}\in X_{h}$ to be the solution of the
following variational
formulation : $\forall v\in X_{h}$, $u^{m+1}$ satisfies

$\frac{1}{\rho}(\nabla(u^{m+1}-u^{m}),\nabla v)+\nu(\nabla u^{m},\nabla
v)+b(u^{m};u^{m+1},v)-(\nabla\cdot v,p^{m})=({f},v),$ (2.7)

and $p^{m+1}\in Q_{h}$ to be the solution of the following variational
formulation : $\forall q\in Q_{h}$, $p^{m+1}$ satisfies

$\alpha((p^{m+1}-p^{m}),q)+\rho(\nabla\cdot u^{m+1},q)=0.$ (2.8)

######  Remark 2.2.

The well-posedness and convergence of the AH scheme was shown in [37], under
some assumptions on the parameter choices for $\rho$ and $\alpha$. Under
similar choices and additional data restrictions beyond the small data
condition, the AH method was shown to be contractive in [7]. While
contractive, the linear convergence rate is rather hard to decipher from the
analysis, which is rather technical (although still an important step
forward). Indeed the rate could be very close to 1, and in the computations
with the AH method in [7], it appears that it often is.

##  3 Convergence analysis of a grad-div stabilized AH method

In this section we consider the following grad-div stabilized AH method, and
will show that it has improved convergence properties over the usual AH
method.

######  Algorithm 3.1.

Let $\rho,\alpha>0$ be parameters, then

  1. 1.

Let $u^{0}\in X_{h}$ and
$p^{0}\in{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{Q}_{h}}$
be the solution of the Stokes problem (2.6).

  2. 2.

For $m\geq 0$, we define $u^{m+1}\in X_{h}$ to be the solution of: $\forall
v\in X_{h}$, $u^{m+1}$ satisfies

$\frac{1}{\rho}(\nabla(u^{m+1}-u^{m}),\nabla v)+\nu(\nabla u^{m},\nabla
v)+b(u^{m};u^{m+1},v)\\\ +\gamma(\nabla\cdot u^{m+1},\nabla\cdot
v)-(\nabla\cdot v,p^{m})=({f},v)$ (3.1)

where the grad-div parameter $\gamma>0$ is a user selected parameter, and
$p^{m+1}\in Q_{h}$ to be the solution of: $\forall q\in Q_{h}$, $p^{m+1}$
satisfies

$\alpha((p^{m+1}-p^{m}),q)+\rho(\nabla\cdot u^{m+1},q)=0.$ (3.2)

We show in this section how grad-div stabilization can provide improved
convergence for AH . First we show linear convergence through a connection to
the classical iterated Picard penalty method, and then we show the classical
convergence analysis with the grad-div term included. Throughout this section
$u$ and $p$ represent the solution of (2.4), and we assume the small data
condition $\kappa<1$ holds.

###  3.1 Linear convergence of the grad-div stabilized AH iteration via a
connection to the iterated Picard penalty method

In this section we show that with certain discretizations, linear convergence
for the grad-div stabilized AH iteration can be established under a small data
condition. In particular, we consider the case of velocity and pressure spaces
satisfying both the inf-sup stability condition and $\nabla\cdot X_{h}=Q_{h}$.
For example, $(P_{2},P_{1}^{disc})$ Scott-Vogelius elements on Alfeld splits
satisfy this property [22]. We prove that with the right parameters, linear
convergence that is equivalent to that of the iterated  Picard penalty (IPP)
method is achieved (which in practice is close to that of the well known
Picard method [14, 9, 27, 33]).

The IPP method for the steady NSE is a classical method that has been well
studied and extensively used [17, 27, 35, 9], and is defined [9] by: Given
$u_{k}\in X_{h},\ p_{k}\in Q_{h}$, solve for $u_{k+1}\in X_{h},\ p_{k+1}\in
Q_{h}$ satisfying

$\displaystyle
b^{*}(u^{m},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}u^{m+1}},v)-(p^{m+1},\nabla\cdot
v)+\nu(\nabla u^{m+1},\nabla v)$ $\displaystyle=(f,v)\ \forall v\in X_{h},$
(3.3) $\displaystyle\epsilon(p^{m+1},q)+(\nabla\cdot u^{m+1},q)$
$\displaystyle=\epsilon(p^{m},q)\ \forall q\in Q_{h},$ (3.4)

where $\epsilon>0$ is a penalty parameter which is generally taken small. It
is proven by Codina in [9] that if the penalty parameter
$\epsilon<\frac{\nu\beta^{2}}{M^{2}}$ and  small data condition $\kappa<1$
holds then both $\|p^{m+1}-p\|$ and $\|\nabla(u^{m}-u)\|$ converge linearly to
$0$, with rate at most

$\mbox{rate}_{IPP}=\left(\frac{1}{2}+\frac{1}{2}\left(1+\frac{2}{r}\right)^{1/2}\right)\left(\kappa+\epsilon\nu^{-1}\beta^{-2}(MC_{1}+r\nu)^{2}\right),$
(3.5)

where

$C_{1}=\left(\nu^{-1}\|f\|_{-1}+\epsilon^{1/2}\nu^{-1/2}\left(\alpha\nu\beta^{-1}\|\nabla(u-u^{0})\|+\|p-p^{0}\|+\|p\|\right)\right)/\left(1-\epsilon^{1/2}\nu^{-1/2}\beta^{-1}M\right)$

and
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}r\geq
2}$ (the optimal $r\geq 2$ will depend on the other parameters). While this is
a complex expression, a rate closer to $\kappa$ –the convergence rate of the
usual Picard iteration [14] which is recovered when $\epsilon=0$– is typically
observed. However, even with larger penalty parameter such as $\epsilon=1$,
IPP enhanced with Anderson acceleration can be very effective, even for larger
data [33].

We now establish that if $\gamma=\rho\alpha^{-1}$, $\rho=\nu^{-1}$ and
$\alpha=\frac{\epsilon}{\nu}$, then the grad-div stabilized AH method is
identical to IPP. From (3.2), since $\nabla\cdot X_{h}=Q_{h}$, we observe that
$p^{m+1}\equiv p^{m}-\rho\alpha^{-1}\nabla\cdot u^{m+1}\equiv
p^{m}-\gamma\nabla\cdot u^{m+1}$, by taking $\gamma=\rho\alpha^{-1}$.
Substituting into (3.1), we obtain

$\frac{1}{\rho}(\nabla(u^{m+1}-u^{m}),\nabla v)+\nu(\nabla u^{m},\nabla
v)+b^{*}(u^{m},u^{m+1},v)-(p^{m+1},\nabla\cdot v)=(f,v),$

and so setting $\rho=\nu^{-1}$ yields

$\nu(\nabla u^{m+1},\nabla v)+b^{*}(u^{m},u^{m+1},v)-(p^{m+1},\nabla\cdot
v)=(f,v).$

Note that the grad-div stabilized AH momentum equation now exactly matches
with the IPP momentum equation (3.3). Substituting $\rho=\nu^{-1}$ into (3.2)
we obtain

$\nu\alpha(p^{m+1}-p^{m},q)+(\nabla\cdot u^{m+1},q)=0.$

Finally setting $\alpha=\frac{\epsilon}{\nu}$, we recover (3.4), thus
establishment of the grad-div stabilized AH and IPP methods are equivalent
when parameters are chosen so that $\gamma=\rho\alpha^{-1}$, $\rho=\nu^{-1}$
and $\alpha=\frac{\epsilon}{\nu}$. With this connection, we have proved the
following theorem for grad-div stabilized AH method.

######  Theorem 3.2.

Suppose grad-div stabilized AH is computed with parameters $\rho=\nu^{-1}$,
$\alpha=\frac{\epsilon}{\nu}$ and $\gamma=\rho\alpha^{-1}=\epsilon^{-1}$, with
user selected penalty parameter $\epsilon<\frac{\nu\beta^{2}}{M^{2}}$. Then,
under the small data condition $\kappa<1$, grad-div stabilized AH method
converges linearly with a rate at most rateIPP defined in (3.5).

######  Remark 3.3.

Just as with IPP, the convergence rate of $\kappa$ will typically be observed
in practice for sufficiently small $\epsilon$. Moreover, with AA, larger
penalty parameters such as $\epsilon=1$ can be used and yield a very efficient
and effective iteration even for $\kappa>1$ as shown in [33].

######  Remark 3.4.

While we do not prove it in our current manuscript, we expect that the
parameters ‘near’ those in the theorem will still provide a linear convergence
by continuity. We believe that it can be proved by using the similar argument
followed in [9], however the additional terms will create a quite challenging
theory  that is beyond the scope of this paper. Our numerical tests show that
grad-div stabilized AH is effective with a rather wide range of parameter
choices, especially for smaller $\epsilon$ and with AA.

###  3.2 Improvement to classical analysis of AH method via grad-div
stabilization

We now consider the improvements to the classical convergence arguments for AH
method, without assuming certain choices of finite elements or meshes other
than $X_{h}\times Q_{h}$ which satisfy the inf-sup condition. Begin the
analysis by adding and subtracting the true solution $(u,p)$ from (3.1) -
(3.2), and then subtracting (2.4) yields

$\frac{1}{\rho}(\nabla(u^{m+1}-u+u-u^{m}),\nabla v)+\nu(\nabla(u^{m}-u),\nabla
v)+b^{*}(u^{m},u^{m+1},v)\\\
-b^{*}(u,u,v)+\gamma(\nabla\cdot(u^{m+1}-u),\nabla\cdot v)-(\nabla\cdot
v,p^{m}-p)=0,$ (3.6)

and

$\alpha((p^{m+1}-p+p-p^{m}),q)+\rho(\nabla\cdot(u^{m+1}-u+u),q)=0,$ (3.7)

respectively.

Denote $e^{m+1}=u^{m+1}-u$ and $e_{p}^{m+1}=p^{m+1}-p$. Taking $v=e^{m+1}$ and
$q=e_{p}^{m+1}$ in (3.6)-(3.7) respectively gives us

$\frac{1}{\rho}(\nabla(e^{m+1}-e^{m}),\nabla e^{m+1})+\nu(\nabla e^{m},\nabla
e^{m+1})+b^{*}(u^{m},u^{m+1},e^{m+1})\\\
-b^{*}(u,u,e^{m+1})+\gamma(\nabla\cdot e^{m+1},\nabla\cdot
e^{m+1})-(\nabla\cdot e^{m+1},e_{p}^{m})=0,$ (3.8)
$\alpha((e_{p}^{m+1}-e_{p}^{m}),e_{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}p}^{m+1})+\rho(\nabla\cdot
e^{m+1},e_{p}^{m+1})=0.$ (3.9)

Since $b^{*}$ is skew-symmetric with respect to its last two arguments,

$\displaystyle b^{*}(u^{m},u^{m+1},e^{m+1})-b^{*}(u,u,e^{m+1})$
$\displaystyle=b^{*}(u^{m},u^{m+1},u^{m+1})-b^{*}(u^{m},u^{m+1},u)-b^{*}(u,u,u^{m+1}),$
$\displaystyle=b^{*}(e^{m},u,u^{m+1}),$
$\displaystyle=b^{*}(e^{m},u,u^{m+1})-b^{*}(e^{m},u,u),$
$\displaystyle=b^{*}(e^{m},u,e^{m+1}).$

Next adding (3.8)-(3.9), using the polarization identity, rearranging and
simplifying the terms, we get

$\frac{1}{2\rho}(\|\nabla(e^{m+1}-e^{m})\|^{2}+\|\nabla
e^{m+1}\|^{2})+\nu\|\nabla e^{m+1}\|^{2}+\gamma\|\nabla\cdot
e^{m+1}\|^{2}+\frac{\alpha}{2\rho}(\|e_{p}^{m+1}-e_{p}^{m}\|^{2}+\|e_{p}^{m+1}\|^{2})\\\
=\frac{1}{2\rho}\|\nabla
e^{m}\|^{2}+\frac{\alpha}{2\rho}\|e_{p}^{m}\|^{2}+\nu(\nabla(e^{m+1}-e^{m}),\nabla
e^{m+1})-b^{*}(e^{m},u,e^{m+1})+(\nabla\cdot e^{m+1},e_{p}^{m}-e_{p}^{m+1}).$
(3.10)

To bound the nonlinear term on the right hand side of last relation, we use
(2.3) and (2.5):

$|-b^{*}(e^{m},u,e^{m+1})|\leq M\|\nabla e^{m}\|\|\nabla u\|\|\nabla
e^{m+1}\|\leq M\nu^{-1}\|f\|_{-1}\|\nabla e^{m}\|\|\nabla
e^{m+1}\|=\kappa\nu\|\nabla e^{m}\|\|\nabla e^{m+1}\|.$

For the third and the last terms in (3.10), we utilize Cauchy-Schwarz and
Young’s inequalities via

$\nu(\nabla(e^{m+1}-e^{m}),\nabla e^{m+1})+(\nabla\cdot
e^{m+1},e_{p}^{m}-e_{p}^{m+1})\\\
\leq\frac{1}{4\rho}\|\nabla(e^{m+1}-e^{m})\|^{2}+\rho\nu^{2}\|\nabla
e^{m+1}\|^{2}+\frac{\alpha}{4\rho}\|e_{p}^{m}-e_{p}^{m+1}\|^{2}+\frac{\rho}{\alpha}\|\nabla\cdot
e^{m+1}\|^{2}.$ (3.11)

Using Young’s and the triangle inequalities provides

$\kappa\nu\|\nabla e^{m}\|\|\nabla
e^{m+1}\|\leq\kappa\nu(\|\nabla(e^{m+1}-e^{m})\|+\|\nabla e^{m+1}\|)\|\nabla
e^{m+1}\|\leq\frac{\kappa\nu}{2}\|\nabla(e^{m+1}-e^{m})\|^{2}+\frac{3\kappa\nu}{2}\|\nabla
e^{m+1}\|^{2},$

and now combining the above bounds, we obtain

$\left(\frac{1}{2\rho}-\frac{\kappa\nu}{2}\right)\|\nabla(e^{m+1}-e^{m})\|^{2}+\left(\frac{1}{2\rho}+\nu-\rho\nu^{2}-\frac{3\kappa\nu}{2}\right)\|\nabla
e^{m+1}\|^{2}\\\ +\left(\gamma-\frac{\rho}{\alpha}\right)\|\nabla\cdot
e^{m+1}\|^{2}+\frac{\alpha}{2\rho}\|e_{p}^{m+1}\|^{2}+\frac{\alpha}{4\rho}\|e_{p}^{m}-e_{p}^{m+1}\|^{2}\leq\frac{1}{2\rho}\|\nabla
e^{m}\|^{2}+\frac{\alpha}{2\rho}\|e_{p}^{m}\|^{2}.$ (3.12)

Provided that $\rho\leq\nu^{-1}\max\left\\{\kappa^{-1},1\right\\}$ and
$\gamma\geq\frac{\rho}{\alpha}$ along with the additional small data
assumption $\kappa<\frac{2}{3}$, the estimate (3.12) is sufficient to provide
convergence of the grad-div stabilized AH method. Comparing to analysis
without the grad-div term from [7], we observe that with grad-div the
coefficient of the left hand side term $\|\nabla e^{m+1}\|^{2}$ is larger and
there are less restrictions on the parameters (including no restriction now on
$\alpha$). The key difference arises from utilizing the left hand side term
$\|\nabla\cdot e^{m+1}\|$, allowing for a larger coefficient of $\|\nabla\cdot
e^{m+1}\|$. Since this is not a proof of contraction, it offers less of a
comparison of rates than the previous section did. Note that if $\nabla\cdot
X_{h}\subset Q_{h}$ then we can follow the proof of [7] to prove a
contraction, however these results would be similar to that of [7] and not
nearly as strong as what is proven above in Theorem 3.2.

##  4 Anderson Acceleration applied to the grad-div stabilized AH Method

In this section, we show that the Anderson acceleration (AA) method can be
applied to the grad-div stabilized AH method (Algorithm 3.1) and will improve
its linear convergence rate. We begin this section with a review of AA and
recent theoretical results. We will proceed to show how the grad-div
stabilized AH method fits into this framework, which in turn allows for
invoking the AA theory. Throughout this section, we assume that the data is
sufficiently small and parameters are chosen so that grad-div stabilized AH
provides a contractive iteration; we specify this assumption below precisely,
after we give some notation.

###  4.1 Anderson acceleration

In this subsection, we provide AA procedure and its convergence properties.
Consider a fixed-point operator $g:Y\rightarrow Y$ where Y is a Hilbert space
equipped with induced norm $\|\cdot\|_{Y}$, and denote
$w_{j}=g(x_{j-1})-x_{j-1}$ as the nonlinear residual, also sometimes is called
the update step. Then, the AA algorithm with depth $m$ (if $m=0$, it returns
to usual Picard iteration) applied to the fixed-point problem $g(x)=x$, reads
as follows.

######  Algorithm 4.1.

(Anderson acceleration with depth $m$ and damping factors $\beta_{k}$)
Step 0: Choose $x_{0}\in Y.$
Step 1: Find $w_{1}\in Y$ such that $w_{1}=g(x_{0})-x_{0}$. Set
$x_{1}=x_{0}+w_{1}$.
Step $k$: For $k=2,3,\ldots$ Set $m_{k}=\min\\{k-1,m\\}.$
[a.] Find $w_{k}=g(x_{k-1})-x_{k-1}$.
[b.] Solve the minimization problem for the Anderson coefficients
$\\{\alpha_{j}^{k}\\}_{k-m_{k}}^{k-1}$

$\displaystyle\\{\alpha_{j}^{k}\\}_{k-m_{k}}^{k-1}=\textstyle\text{argmin}\left\|\left(1-\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k}\right)w_{k}+\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k}w_{k-j}\right\|_{Y}.$
(4.1)

[c.] For damping factor $0<\beta_{k}\leq 1$, set

$\displaystyle\textstyle
x_{k}=(1-\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k})x_{k-1}+\sum_{j=k-m_{k}}^{k-1}\alpha_{j}^{k}x_{j-1}+\beta_{k}\left((1-\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k})w_{k}+\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k}w_{k-j}\right).$
(4.2)

To understand how AA improves convergence, we define the optimization gain
factor $\theta_{k}$ by

$\displaystyle\theta_{k}=\frac{\left\|{(1-\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k})w_{k}+\sum\limits_{j=k-m_{k}}^{k-1}\alpha_{j}^{k}w_{k-j}}\right\|_{Y}}{\|w_{k}\|_{Y}},$

which characterizes the improvement in fixed-point convergence rate as
proposed in [30, 31].

The following assumptions from [30] provide sufficient conditions on the fixed
point operator $g$ for the convergence and acceleration results.

######  Assumption 4.2.

Assume $g\in C^{1}(Y)$ has a fixed point $x^{\ast}$ in $Y$, and there are
positive constants $C_{0}$ and $C_{1}$ with

  1. 1.

$\|g^{\prime}(x)\|_{Y}\leq C_{0}$ for all $x\in Y$, and

  2. 2.

$\|g^{\prime}(x)-g^{\prime}(y)\|_{Y}\leq C_{1}\|x-y\|_{Y}$ for all $x,y\in Y$.

######  Assumption 4.3.

Assume there is a constant $\sigma>0$ for which the differences between
consecutive residuals and iterates satisfy

$\displaystyle\|w_{{k}+1}-w_{k}\|_{Y}\geq\sigma\|x_{k}-x_{{k}-1}\|_{Y},\quad{k}\geq
1.$ (4.3)

Assumption 4.2 will be verified for Picard fixed-point operator grad-div
stabilized AH method in next sections. Also, Assumption 4.3 can be verified
easily for this method which is contractive under small data and particular
parameter choices. Under Assumptions 4.2 and 4.3, the following result from
[30], generates a bound on the residual $\|w_{k+1}\|$ in terms of the previous
residual $\|w_{k}\|$.

######  Theorem 4.4 (Pollock et al., 2021).

Let Assumptions 4.2 and 4.3 hold, and suppose the direction sines between each
column $j$ of matrix

$\displaystyle F_{j}$ $\displaystyle:=\left((w_{j}-w_{j-1})(w_{j-1}-w_{j-2})\
\dotsc\ (w_{j-m_{j}+1}-w_{j-m_{j}})\right)$

and the subspace spanned by the preceeding columns satisfy
$|\sin(f_{j,i},\text{span }\\{f_{j,1},\ldots,f_{j,i-1}\\})|\geq c_{s}>0$, for
$j=k-m_{k},\ldots,k-1$. Then the residual $w_{k+1}=g(x_{k})-x_{k}$ from
Algorithm 4.1 (depth $m$) satisfies the following bound.

$\displaystyle\left\|{w_{k+1}}\right\|$
$\displaystyle\leq\left\|{w_{k}}\right\|\Bigg{(}\theta_{k}((1-\beta_{k})+C_{0}\beta_{k})+\frac{CC_{1}\sqrt{1-\theta_{k}^{2}}}{2}\bigg{(}\left\|{w_{k}}\right\|h(\theta_{k})$
$\displaystyle+2\sum_{n=k-{m_{k}}+1}^{k-1}(k-n)\left\|{w_{n}}\right\|h(\theta_{n})+m_{k}\left\|{w_{k-m_{k}}}\right\|h(\theta_{k-m_{k}})\bigg{)}\Bigg{)},$
(4.4)

where each $h(\theta_{j})\leq C\sqrt{1-\theta_{j}^{2}}+\beta_{j}\theta_{j}$,
and $C$ depends on $c_{s}$ and the implied upper bound on the direction
cosines.

In (4.4), the optimization gain $\theta_{k}$ is scaling the first-order term,
which is residual in the standard fixed-point iteration. On the other hand,
the higher-order terms are scaled by a factor of $\sqrt{1-\theta_{k}^{2}}$,
which implies if the optimization works, the relative weight of the higher-
order terms increase, otherwise the relative weight of the first-order term
increase in (4.4).

###  4.2 Grad-div stabilized AH method as a fixed point iteration

In this subsection, we define the fixed-point operator $G$ which is associated
with grad-div stabilized AH iteration. Note that in this section, $u,p$ are
generic functions and are not the steady NSE solution as in the previous
section.

######  Definition 4.5.

Define mapping $G:(X_{h},Q_{h})\rightarrow(X_{h},Q_{h})$,
$G(u,p)=(G_{1}(u,p),G_{2}(u,p))$ such that for any $(v,q)\in(X_{h},Q_{h})$

$\displaystyle\rho^{-1}(\nabla(G_{1}(u,p)-u),\nabla v)+\nu(\nabla u,\nabla
v)+b(u;G_{1}(u,p),v)$ $\displaystyle+\gamma(\nabla\cdot G_{1}(u,p),\nabla\cdot
v)-(p,\nabla\cdot v)$ $\displaystyle=(f,v),$ (4.5)
$\displaystyle\alpha(G_{2}(u,p)-p,q)+\rho(\nabla\cdot G_{1}(u,p),q)$
$\displaystyle=0.$ (4.6)

Now, we show that $G$ is well-defined and bounded with respect to the norm
$\|(v,q)\|_{H}\eqqcolon\sqrt{\|\nabla v\|^{2}+\alpha\|q\|^{2}}$ on
$X_{h}\times Q_{h}$.

######  Lemma 4.6.

The operator $G$ is well-defined. Moreover, if
$G(u,p)=(G_{1}(u,p),G_{2}(u,p))$ is the solution of (4.5)-(4.6), then the
following inequality holds

$\displaystyle\|G(u,p)\|_{H}\leq\sqrt{\frac{(1+4\rho^{2}\nu^{2})}{N}}\|\nabla
u\|+\sqrt{\frac{\alpha}{N}}\|p\|+\sqrt{\frac{4\rho^{2}}{N}}\|f\|_{-1},$ (4.7)

where $N\coloneqq\min\\{\frac{1}{2}-\alpha^{-1}\rho^{2},1\\}$.

###### Proof.

Assume that a solution exists. Then, choosing $v=G_{1}(u,p)$ and
$q=G_{2}(u,p)$ eliminates the nonlinear term and yields

$\displaystyle\frac{1}{2}\left(\|\nabla G_{1}(u,p)\|^{2}-\|\nabla
u\|^{2}+\|\nabla(G_{1}(u,p)-u)\|^{2}\right)+\rho\gamma\|\nabla\cdot
G_{1}(u,p)\|^{2}-\rho(p,\nabla\cdot G_{1}(u,p))$
$\displaystyle=\rho(f,G_{1}(u,p))-\rho\nu(\nabla u,\nabla G_{1}(u,p)),$
$\displaystyle\frac{\alpha}{2}\left(\|G_{2}(u,p)\|^{2}-\|p\|^{2}+\|G_{2}(u,p)-p\|^{2}\right)+\rho(\nabla\cdot
G_{1}(u,p),G_{2}(u,p))=0.$

thanks to the polarization identity. Then, combining the above equations gives
that

$\displaystyle\frac{1}{2}\left(\|\nabla
G_{1}(u,p)\|^{2}+\|\nabla(G_{1}(u,p)-u)\|^{2}\right)+\frac{\alpha}{2}\left(\|G_{2}(u,p)\|^{2}+\|G_{2}(u,p)-p\|^{2}\right)+\rho\gamma\|\nabla\cdot
G_{1}(u,p)\|^{2}$ $\displaystyle=\frac{1}{2}\|\nabla
u\|^{2}+\frac{\alpha}{2}\|p\|^{2}+\rho(f,G_{1}(u,p))-\rho\nu(\nabla u,\nabla
G_{1}(u,p))-\rho(G_{2}(u,p)-p,\nabla\cdot G_{1}(u,p)).$

Multiplying the both sides of the last relation by $2$, dropping positive
terms $\|\nabla(G_{1}(u,p)-u)\|^{2}$, $\rho\gamma\|\nabla\cdot
G_{1}(u,p)\|^{2}$, and $\|G_{2}(u,p)-p\|^{2}$ on the left hand side, and using
Hölder’s and Young’s inequalities produce

$\displaystyle\left(\frac{1}{2}-\alpha^{-1}\rho^{2}\right)\|\nabla
G_{1}(u,p)\|^{2}+\alpha\|G_{2}(u,p)\|^{2}\leq(1+4\rho^{2}\nu^{2})\|\nabla
u\|^{2}+\alpha\|p\|^{2}+\rho^{2}4\|f\|^{2}_{-1}.$

Letting $N\coloneqq\min\\{\frac{1}{2}-\alpha^{-1}\rho^{2},1\\}$, and dividing
both sides by $N$, we get

$\displaystyle\|G(u,p)\|^{2}\leq\frac{(1+4\rho^{2}\nu^{2})}{N}\|\nabla
u\|^{2}+\frac{\alpha}{N}\|p\|^{2}+\frac{4\rho^{2}}{N}\|f\|^{2}_{-1}.$

Then, taking the square root of both sides reduces it to (4.7). Since $G$ is
linear and finite dimensional, showing that the solution $G(u,p)$ is bounded
continuously by the data implies solution uniqueness and thus existence as
well.

∎

We now rewrite the grad-div stabilized AH method in terms of a mapping
$G:(X_{h},Q_{h})\rightarrow(X_{h},Q_{h})$ that satisfies for $m\geq 0$

$G(u^{m},p^{m})=(G_{1}(u^{m},p^{m}),G_{2}(u^{m},p^{m})):=(u^{m+1},p^{m+1}),$

where $(u^{m},p^{m})$ is the $m^{th}$ iteration of the A-H method described in
Algorithm 3.1.

###  4.3 Applying AA to the grad-div stabilized AH iteration

In this subsection, we show the sufficient smoothness properties of the
associated fixed point operator $G$ for the grad-div stabilized AH iteration
to apply AA theory. Now, we show Lipschitz continuity of $G$.

######  Lemma 4.7.

For any $(u,p),(w,z)\in(X_{h},Q_{h})$, we have

$\displaystyle\|G(u,p)-G(w,z)\|_{H}\leq C_{L}\|(u,p)-(w,z)\|_{H}$ (4.8)

where
$C_{L}=\max\\{\sqrt{\frac{2}{K}},\left(\frac{1-\rho\nu}{\sqrt{K}}+\frac{\rho}{\sqrt{KN}}\left((1+2\rho\nu)\|\nabla
u\|+\sqrt{\alpha}\|p\|+2\rho\|f\|_{-1}\right)\right)\\}$.

######  Remark 4.8.

Note that we have already discussed that in case of sufficiently small data
and particularly chosen parameters, Algorithm 3.1 is contractive. In this
section we assume to be in the contractive setting, and thus we have that the
Lipschitz constant $C_{L}$ in Lemma 4.7 is less than $1$.

###### Proof.

Subtracting (4.5) with $(w,z)$ from (4.5) with $(u,p)$ gives

$\displaystyle(\nabla(G_{1}(u,p)-G_{1}(w,z)),\nabla v)-(\nabla(u-w),\nabla
v)+\rho\nu(\nabla(u-w),\nabla v)+\rho b(u;$ $\displaystyle
G_{1}(u,p)-G_{1}(w,z),v)$ $\displaystyle+\rho
b(u-w;G_{1}(w,z),v)+\rho\gamma(\nabla\cdot(G_{1}(u,p)-G_{1}(w,z)),\nabla\cdot
v)$ $\displaystyle=\rho(p-z,\nabla\cdot v),$
$\displaystyle\alpha(G_{2}(u,p)-G_{2}(w,z),q)-\alpha(p-z,q)+\rho(\nabla\cdot(G_{1}(u,p)-G_{1}(w,z)),q)$
$\displaystyle=0.$

Then, setting $v=G_{1}(u,p)-G_{1}(w,z)$ and $q=G_{2}(u,p)-G_{2}(w,z)$ which
eliminates the first nonlinear term on the left hand side of the first
equation, and combining these equations provide

$\displaystyle\|\nabla(G_{1}(u,p)-G_{1}(w,z))\|^{2}+\alpha\|G_{2}(u,p)-G_{2}(w,z)\|^{2}+\rho\gamma\|\nabla\cdot(G_{1}(u,p)-G_{1}(w,z))\|^{2}$
$\displaystyle=\rho(p-z,\nabla\cdot(G_{1}(u,p)-G_{1}(w,z)))+(1-\rho\nu)(\nabla(u-w),\nabla(G_{1}(u,p)-G_{1}(w,z)))$
$\displaystyle-\rho
b(u-w;G_{1}(w,z),G_{1}(u,p)-G_{1}(w,z))+\alpha(p-z,G_{2}(u,p)-G_{2}(w,z))$
$\displaystyle-\rho(\nabla\cdot(G_{1}(u,p)-G_{1}(w,z)),G_{2}(u,p)-G_{2}(w,z)).$

By dropping positive term
$\rho\gamma\|\nabla\cdot(G_{1}(u,p)-G_{1}(w,z))\|^{2}$ on the left hand side
and using Cauchy-Schwarz and Young’s inequalities, Lemma 4.6 and (2.3) , we
get

$\displaystyle\left(\frac{1}{2}-\frac{3\rho^{2}\alpha^{-1}}{4}\right)\|\nabla(G_{1}(u,p)-G_{1}(w,z))\|^{2}+\frac{\alpha}{4}\|G_{2}(u,p)-G_{2}(w,z)\|^{2}\leq
2\alpha\|p-z\|^{2}$
$\displaystyle+\left((1-\rho\nu)^{2}+\rho^{2}\left(\frac{(1+4\rho^{2}\nu^{2})}{N}\|\nabla
u\|^{2}+\frac{\alpha}{N}\|p\|^{2}+\frac{4\rho^{2}}{N}\|f\|^{2}_{-1}\right)\right)\|\nabla(u-w)\|^{2}.$

Defining
$K=\min\\{\left(\frac{1}{2}-\frac{3\rho^{2}\alpha^{-1}}{4}\right),\frac{1}{4}\\}$
and then dividing both sides by $K$ yields (4.8). ∎

Next, we define an operator $G^{\prime}$ and show it is indeed the Fréchet
derivative of the operator of $G.$

######  Definition 4.9.

Given $(u,p)\in X_{h}\times Q_{h}$, define an operator
$G^{\prime}(u,p;\cdot,\cdot):X_{h}\times Q_{h}\rightarrow X_{h}\times Q_{h}$
by

$\displaystyle
G^{\prime}(u,p;w,s)\eqqcolon(G^{\prime}_{1}(u,p;w,s),G^{\prime}_{2}(u,p;w,s))$

satisfying for all $(w,s)\in X_{h}\times Q_{h}$.

$\displaystyle(\nabla(G^{\prime}_{1}(u,p;w,s)-w),\nabla v)+\rho\nu(\nabla
w,\nabla v)+$ $\displaystyle\rho b(w;G_{1}(u,p),v)+\rho
b(u;G^{\prime}_{1}(u,p;w,s),v)$ (4.9) $\displaystyle\rho\gamma(\nabla\cdot
G^{\prime}_{1}(u,p;w,s),\nabla\cdot v)-\rho(s,\nabla\cdot v)$
$\displaystyle=0,$
$\displaystyle\alpha(G^{\prime}_{2}(u,p;w,s)-s,q)+\rho(\nabla\cdot
G^{\prime}_{1}(u,p;w,s),q)$ $\displaystyle=0.$

######  Lemma 4.10.

The operator $G^{\prime}$ is well-defined for all $(u,p),(w,s)\in X_{h}\times
Q_{h}$ such that

$\displaystyle\|G^{\prime}(u,p;w,s)\|_{H}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}C_{L}}\|(w,s)\|_{H}$
(4.10)

###### Proof.

Adding equations in (4.9) and setting $v=G^{\prime}_{1}(u,p;w,s)$ and
$q=G^{\prime}_{2}(u,p;w,s)$ produces

$\displaystyle\|\nabla G^{\prime}_{1}(u,p;w,s)\|^{2}+\rho\gamma\|\nabla\cdot
G^{\prime}_{1}(u,p;w,s)\|^{2}+\alpha\|G^{\prime}_{2}(u,p;w,s)\|^{2}$ (4.11)
$\displaystyle=(1-\rho\nu)(\nabla w,\nabla
G^{\prime}_{1}(u,p;w,s))+\alpha(s,G^{\prime}_{2}(u,p;w,s))-\rho
b(w;G_{1}(u,p),G^{\prime}_{1}(u,p;w,s))$ $\displaystyle+\rho(s,\nabla\cdot
G^{\prime}_{1}(u,p;w,s))-\rho(\nabla\cdot
G^{\prime}_{1}(u,p;w,s),G^{\prime}_{2}(u,p;w,s)).$

Then, by dropping positive term $\rho\gamma\|\nabla\cdot
G^{\prime}_{1}(u,p;w,s)\|^{2}$ on the left hand side, applying Cauchy-Schwarz
and Young’s inequalities produces

$\displaystyle\|\nabla
G^{\prime}_{1}(u,p;h,s)\|^{2}+\alpha\|G^{\prime}_{2}(u,p;h,s)\|^{2}$
$\displaystyle\leq(1-\rho\nu)^{2}\|\nabla h\|^{2}+\frac{1}{4}\|\nabla
G^{\prime}_{1}(u,p;h,s)\|^{2}+\alpha\|s\|^{2}+\frac{\alpha}{4}\|G^{\prime}_{2}(u,p;h,s)\|^{2}$
$\displaystyle+\rho^{2}M^{2}\left(\frac{(1+4\rho^{2}\nu^{2})}{N}\|\nabla
u\|^{2}+\frac{\alpha}{N}\|p\|^{2}+\frac{4\rho^{2}}{N}\|f\|^{2}_{-1}\right)\|\nabla
h\|^{2}+\frac{1}{4}\|\nabla G^{\prime}_{1}(u,p;h,s)\|^{2}$
$\displaystyle+\alpha\|s\|^{2}+\frac{\rho^{2}\alpha^{-1}}{4}\|\nabla
G^{\prime}_{1}(u,p;h,s)\|^{2}+\frac{\rho^{2}\alpha^{-1}}{2}\|\nabla
G^{\prime}_{1}(u,p;h,s)|^{2}+\frac{\alpha}{2}\|G^{\prime}_{2}(u,p;h,s)\|^{2}$

thanks to the Lemma 4.6 and (2.3). After rearranging the terms, we get

$\displaystyle\left(\frac{1}{2}-\frac{3\rho^{2}\alpha^{-1}}{4}\right)\|\nabla
G^{\prime}_{1}(u,p;h,s)\|^{2}+\frac{\alpha}{4}\|G^{\prime}_{2}(u,p;h,s)\|^{2}$
$\displaystyle\leq\left((1-\rho\nu)^{2}+\rho^{2}M^{2}\left(\frac{(1+4\rho^{2}\nu^{2})}{N}\|\nabla
u\|^{2}+\frac{\alpha}{N}\|p\|^{2}+\frac{4\rho^{2}}{N}\|f\|^{2}_{-1}\right)\right)\|\nabla
h\|^{2}+2\alpha\|s\|^{2}$

Taking the square root of both sides and considering constant $K$ which is
defined in proof of Lemma 4.7 finishes the proof.

∎

Now, we prove that $G^{\prime}$ is Fréchet derivative operator of $G.$

######  Lemma 4.11.

For any $(u,p)\in X_{h}\times Q_{h}$

$\displaystyle\|\ G(u+w,p+s)-G(u,p)-G^{\prime}(u,p;w,s)\|_{H}\leq\frac{\rho
MC_{L}}{\sqrt{1-\rho\alpha^{-1}}}\|(w,s)\|_{H}^{2}.$ (4.12)

###### Proof.

Denote
$\eta_{1}=G_{1}(u+w,p+s)-G_{1}(u,p)-G_{1}^{\prime}(u,p;w,s),\eta_{2}=G_{2}(u+w,p+s)-G_{2}(u,p)-G_{2}^{\prime}(u,p;w,s).$
Subtracting the sum of (4.9) and (4.5) from the equation (4.5) with
$(u+w,p+s)$ yields

$\displaystyle(\nabla\eta_{1},\nabla
v)+\rho\gamma(\nabla\cdot\eta_{1},\nabla\cdot v)+\rho b(u;\eta_{1},v)+\rho
b(w;G_{1}(u+w,p+s)-G_{1}(u,p),v)$ $\displaystyle=0,$
$\displaystyle\alpha(\eta_{2},q)+\rho(\nabla\cdot\eta_{1},q)$
$\displaystyle=0.$

Choosing $v=\eta_{1}$ and $q=\eta_{2}$ vanishes first nonlinear term on the
left hand side, and combining these equations gives

$\|\nabla\eta_{1}\|^{2}+\alpha\|\eta_{2}\|^{2}+\rho\gamma\|\nabla\cdot\eta_{1}\|^{2}=-\rho(\nabla\cdot\eta_{1},\eta_{2})-\rho
b(w;G_{1}(u+w,p+s)-G_{1}(u,p),\eta_{1}).$ (4.13)

Applying Cauchy-Schwarz on the right hand side and dropping
$\rho\gamma\|\nabla\cdot\eta_{1}\|^{2}$ on the left hand side yields

$\displaystyle\|\nabla\eta_{1}\|^{2}+\alpha\|\eta_{2}\|^{2}\leq\rho\|\nabla\cdot\eta_{1}\|\|\eta_{2}\|+M\rho\|\nabla
w\|\|\nabla G_{1}(u+w,p+s)-G_{1}(u,p)\|\|\nabla\eta_{1}\|.$

Now with Young’s inequality and Lemma 4.7, we obtain using
$\|\cdot\nabla\|\leq\|\nabla\|$ and $\|\nabla w\|^{2}\leq\|\nabla
w\|^{2}+\alpha\|s\|^{2}=\|(w,s)\|_{H}^{2}$ that

$\displaystyle\|\nabla\eta_{1}\|^{2}+\alpha\|\eta_{2}\|^{2}\leq\frac{\rho^{2}\alpha^{-1}}{2}\|\nabla\eta_{1}\|^{2}+\frac{\alpha}{2}\|\eta_{2}\|^{2}+\frac{\rho^{2}}{2}M^{2}C_{L}^{2}\|(w,s)\|^{4}_{H}+\frac{1}{2}\|\nabla\eta_{1}\|^{2}.$

Then, rearrange and obtain the following

$\displaystyle\left(\frac{1}{2}-\frac{\rho^{2}\alpha^{-1}}{2}\right)\|\nabla\eta_{1}\|^{2}+\frac{\alpha}{2}\|\eta_{2}\|^{2}$
$\displaystyle\leq\frac{\rho^{2}}{2}M^{2}C_{L}^{2}\|\nabla
w\|^{2}\|\nabla(G_{1}(u+w,p+s)-G_{1}(u,p))\|^{2}$
$\displaystyle\leq\frac{\rho^{2}}{2}M^{2}C_{L}^{2}\|(w,s)\|_{H}^{4},$

where we used (2.3). Then, applying the definition of $\eta_{1}$ and
$\eta_{2}$, dividing both sides by
$\left(\frac{1}{2}-\frac{\rho^{2}\alpha^{-1}}{2}\right)$ and taking square
roots give that $G^{\prime}$ is indeed the Fréchet derivative of $G$ which
satisfies (4.12). ∎

We now proceed to show that $G^{\prime}$ is Lipschitz continuous over
$X_{h}\times Q_{h}$.

######  Lemma 4.12.

G is Lipschitz continuously differentiable on $X_{h}\times Q_{h}$, such that
for all $u,w,\theta\in X_{h}$ and $p,s,\xi\in Q_{h}$,

$\displaystyle\|G^{\prime}(u+w,p+s;\theta,\xi)-G^{\prime}(u,p;\theta,\xi)\|_{H}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\left(\frac{4\rho^{2}M^{2}C_{L}^{2}}{1-\rho^{2}\alpha^{-1}}\right)^{1/2}}\|(\theta,\xi)\|_{H}\|(w,s)\|_{H}.$
(4.14)

###### Proof.

Subtracting (4.9) with $G^{\prime}(u,p;\theta,\xi)$ from (4.9) with
$G^{\prime}(u+w,p+s;\theta,\xi)$ and denoting
$e_{1}\eqqcolon
G^{\prime}_{1}(u+w,p+s;\theta,\xi)-G^{\prime}_{1}(u,p;\theta,\xi)$ and
$e_{2}\eqqcolon
G^{\prime}_{2}(u+w,p+s;\theta,\xi)-G^{\prime}_{2}(u,p;\theta,\xi)$ yield

$\displaystyle(\nabla e_{1},\nabla v)+\rho
b(\theta;G_{1}(u+w,p+s)-G_{1}(u,p),v)+\rho b(u;e_{1},v)$ $\displaystyle+\rho
b(w;G^{\prime}_{1}(u+w,p+s;\theta,\xi),v)+\rho\gamma(\nabla\cdot
e_{1},\nabla\cdot v)$ $\displaystyle=0,$
$\displaystyle\alpha(e_{2},q)+\rho(\nabla\cdot e_{1},q)$ $\displaystyle=0.$

Setting $v=e_{1}$ and $q=e_{2}$ vanishes the third term on the left hand side
of the first equality and adding these equations provide

$\displaystyle\|\nabla e_{1}\|^{2}+\alpha\|e_{2}\|^{2}+\rho\gamma\|\nabla\cdot
e_{1}\|^{2}$ $\displaystyle=-\rho(\nabla\cdot e_{1},e_{2})-\rho
b(\theta;G_{1}(u+w,p+s)-G_{1}(u,p),e_{1})-\rho
b(w;G^{\prime}_{1}(u+w,p+s;\theta,\xi),e_{1}).$

Dropping the positive term $\rho\gamma\|\nabla\cdot e_{1}\|^{2}$ on the left
hand side, applying Cauchy-Schwarz and Young’s inequalities and (2.3), using
Lemma 4.7 and 4.10, we get

$\displaystyle\left(\frac{1}{2}-\frac{\rho^{2}\alpha^{-1}}{2}\right)\|\nabla
e_{1}\|^{2}+\frac{\alpha}{2}\|e_{2}\|^{2}$
$\displaystyle\leq\rho^{2}M^{2}C_{L}^{2}\|\nabla\theta\|^{2}\|(w,s)\|^{2}_{H}+\rho^{2}M^{2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}C_{L}^{2}}\|\nabla
w\|^{2}\|(\theta,\xi)\|^{2}$
$\displaystyle\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}2\rho^{2}M^{2}C_{L}^{2}}\|(\theta,\xi)\|_{H}^{2}\|(w,s)\|^{2}_{H}.$

Dividing both sides by
$\left(\frac{1}{2}-\frac{\rho^{2}\alpha^{-1}}{2}\right)$ gives

$\displaystyle\|(e_{1},e_{2})\|_{H}^{2}=\|\nabla
e_{1}\|^{2}+\alpha\|e_{2}\|^{2}\leq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\frac{4\rho^{2}M^{2}C_{L}^{2}}{1-\rho^{2}\alpha^{-1}}}\|(\theta,\xi)\|_{H}^{2}\|(w,s)\|^{2}_{H}.$

Then taking the square roots of both sides gives that $G^{\prime}$ is
Lipschitz continuous, and (4.14) holds. ∎

###  4.4 Convergence of the Anderson Accelerated AH algorithm for steady NSE

In previous subsection, we proved that the solution operator $G$ associated
with grad-div stabilized AH iteration (4.5)-(4.6) satisfies Assumption 4.2
which is the one of sufficient conditions to apply the one-step residual bound
of [30]. Also, Assumption 4.3 is satisfied since $G$ is contractive under
small data condition and certain parameter choices.

Under these assumptions and with Lemmas 4.7, 4.11, 4.12 and Theorem 4.4, we
have established the convergence of (4.5)-(4.6) where $G$ is the solution
operator associated with grad-div stabilized AH iteration.

######  Theorem 4.13.

For any step $k>m$ with $\alpha_{m}^{k}\neq 0$, the following bound holds for
the grad-div stabilized AH iteration (4.5)-(4.6)

$\displaystyle\|(w_{k+1},z_{k+1})\|_{H}\leq$
$\displaystyle\theta_{k}(1-\beta_{k}+\beta_{k}C_{L})\|(w_{k},z_{k})\|_{H}$
$\displaystyle+C\sqrt{1-\theta_{k}^{2}}\|(w_{k},z_{k})\|_{H}\sum\limits_{j=1}^{m}\|(w_{k-j+1},z_{k-j+1})\|_{H},$

for the residual $(w_{k},z_{k})$, where $\theta_{k}$ is the gain from the
optimization problem, $C_{L}$ is the Lipschitz constant of $G$ defined in
Lemma 4.7, and $C$ depending on $\theta_{k},\beta_{k},C_{L}$.

This theorem tells us that (4.5)-(4.6), with a good initial guess, converges
linearly with rate $\theta_{k}(1-\beta_{k}+\beta_{k}C_{L})<1$, which improves
on Algorithm 3.1 due to the scaling $\theta_{k}$ and the damping factor
$\beta_{k}$. In the case $G$ is contractive, i.e. $C_{L}<1$, then the optimal
choice for relaxation is $\beta_{k}=1$.

##  5 Numerical Experiments

In this section, we perform several numerical tests to illustrate the theory
above and to show how the grad-div stabilized, Anderson accelerated AH
algorithm can be an effective and efficient solver for the steady NSE. The
stopping criteria for all of our tests is $\|u_{k}-u_{k-1}\|\leq 10^{-6}$.

| $\gamma=0$
---|---
Taylor-Hood |
Scott-Vogelius |

| $\gamma=1$
---|---
Taylor-Hood |
Scott-Vogelius |

Figure 1:  Shown above is convergence behavior for the AH method applied to
the driven cavity problem with varying $Re$ and $\rho$. The top plots are
without grad-div stabilization and the bottom plots are with it. The first and
third rows are for TH elements, and the second and last rows are SV.

###  5.1 Lid-driven cavity

We first test the AH method for steady NSE on the lid-driven cavity problem.
The domain for the problem is the unit square $\Omega=(0,1)^{2}$ and we impose
Dirichlet boundary conditions by $u|_{y=1}=(1,0)^{T}$ and $u=0$ everywhere
else. We choose the parameter $\alpha=\nu^{-1}$. We first illustrate how grad-
div stabilization improves the AH method, and show the dramatic improvement
offered by $(P_{2},P^{disc}_{1})$ Scott-Vogelius (SV) over $(P_{2},P_{1})$
Taylor-Hood (TH). The second test shows even further dramatic improvement by
incorporating AA.

####  5.1.1 The effect of grad-div stabilization and comparison of Scott-
Vogelius vs Taylor-Hood

Our analysis above suggests that convergence of AH will be improved from using
grad-div stabilization, and also from using SV elements instead of TH since
the connection to the iterated Picard penalty method is only made for SV
elements. Hence we now compare the AH method for both element choices, with
and without grad-div stabilization (using parameter $\gamma=1$). We run the
tests for varying $\rho$ (to try to find a good choice of parameter $\rho$)
and with varying $Re=\nu^{-1}$. For these tests a uniform mesh with $h=1/32$
is used.

Results are shown in Figure 1. We observe the best AH results clearly come
from using SV instead of TH, and using $\gamma=1$ with SV gives by far the
best results. In most cases, TH fails to converge for any $\rho$. Based on
these results, we use SV elements for the rest of the numerical tests in this
paper.

####  5.1.2 Anderson accelerated grad-div stabilized AH method with SV
elements

We now consider the same test problem, using the AH iteration only with SV
elements and $\gamma=1$, but now adding AA. We test AA depths $m=0$ (no
acceleration),$1,5$ and $10$. Figure 2 and 3 show convergence results obtained
by Anderson accelerated grad-div stabilized AH method for $Re=100$ and $1000$,
respectively, for varying $\rho$. As depth increases, the number of iterations
decreases significantly. The fastest convergence is obtained with depths
$m=50$, but there is not much improvement past $m=5$.

We also note that for optimally chosen $\rho$ in this setting (i.e. $\rho=20$
for $Re=100$ and $\rho=50$ for $Re=1000$ based on test above) with SV elements
and grad-div stabilization, there is not much difference in convergence from
AA. However, for slightly non optimal $\rho$, there can be a dramatic
improvement from AA. Since one often does not know optimal $\rho$ a priori,
the expected case in practice is using a non-optimal $\rho$.

Comparing to existing literature, for $Re=100$ driven cavity it is reported in
[7] that 731 iterations of AH were needed to converge to the same tolerance
used herein and with Taylor-Hood elements, $\rho=1.2$ and $\alpha=70$. Our
results with Taylor-Hood elements and no grad-div stabilization were similarly
bad, see figure 1 in row 1; in fact, that they got convergence at all for this
test is rather extraordinary. With SV elements and grad-div stabilization,
figure 1 shows that with $\alpha=100$ and $\rho=20$, convergence is achieved
in 80 iterations. With less optimal parameter choices, AA can still keep the
total iteration count low, see figure 2.

Figure 2: Convergence of Anderson accelerated grad-div stabilized Arrow-
Hurwicz method for $Re=100$, different $\rho$’s, $\gamma=1$ with varying $m$.

Figure 3: Convergence of Anderson accelerated grad-div stabilized Arrow-
Hurwicz method for $Re=1000$, different $\rho$’s, $\gamma=1$ with varying $m$.

####  5.1.3 Driven cavity with $Re$=5,000 and $Re$=10,000

As a final test with the driven cavity, we consider the case of $Re$=5,000 and
$Re$=10,000 with $h$=1/64. This is a difficult problem for nonlinear solvers
[31, 30], and we show now that with the right parameter choices, the AH method
can be effective for this problem. For 5,000, $\rho=100,\ \gamma=1,\ m=100$
were used to obtain convergence in 464 iterations. For 10,000, $\rho=150,\
\gamma=10,\ m=100$ was used to obtain convergence in 217 iterations.
Streamlines for both solutions are shown in Figure 4, and they are in good
agreement with those from [13] even though we use a coarser mesh. We note
that, to date, the highest $Re$ for successful lid driven cavity computations
in the literature is 1000 in [7].

Figure 4: Show above are streamlines for $Re=5,000$ and $10,000$ driven cavity
solutions.

###  5.2 Channel flow past a step

Figure 5: Shown above are the mesh for channel flow past a step (left) and
Re=100 velocity solution (right).

For our last test, we consider 2D channel flow past a step with
$Re=\nu^{-1}=100$. The domain for this problem is a $40\times 10$ rectangular
channel, with a $1\times 1$ ‘step’ placed 5 units into the channel at the
bottom. The triangulation we use is shown in Figure 5, along with the $Re=100$
solution found with our solver (which is consistent with solutions from the
literature [21, 15, 16]). The discretization uses $(P_{2},P_{1}^{disc})$ SV
elements that provided 32,682 velocity degrees of freedom.

First we consider $\gamma=10=\varepsilon^{-1}$, noting that obtaining
convergence with $\gamma=1$ proved very difficult and we were not able to find
a parameter set that gave convergence. With $\gamma=10$, we computed four
parameter sets: $(\rho=50,\ \alpha=\frac{\varepsilon}{\nu})$ , $(\rho=50,\
\alpha=\frac{1}{\nu})$, $(\rho=100,\ \alpha=\frac{\varepsilon}{\nu})$ \- which
is exactly the iterated penalty Picard method, and $(\rho=100,\
\alpha=\frac{1}{\nu})$. Convergence plots for each of these parameter sets and
varying $m$ are shown in figure 6, and we observe that $m=100$ is the best
choice for AA for all cases, and that AH with parameters chosen to match IPP
performs significantly worse than other parameter choices. The choice
$\rho=50$ and $\alpha=\nu^{-1}$ with $m=100$ was very effective. Results
improve with $\gamma=100$ for all parameter sets, see figure 7. Here, again
$AH$ improves on IPP, with $\rho=50$ and $\rho=100$ withe $m=100$ were very
effective parameter choices.

$(\rho=50,\ \alpha=\frac{\varepsilon}{\nu})$ $(\rho=50,\
\alpha=\frac{1}{\nu})$ $(\rho=100,\ \alpha=\frac{\varepsilon}{\nu})$ (IPP)
$(\rho=100,\ \alpha=\frac{1}{\nu})$

Figure 6: Convergence of Anderson accelerated grad-div stabilized AH method
and IPP iteration for $Re=100$ for varying parameters and $\gamma=10$.

$(\rho=50,\ \alpha=\frac{\varepsilon}{\nu})$ $(\rho=50,\
\alpha=\frac{1}{\nu})$ $(\rho=100,\ \alpha=\frac{\varepsilon}{\nu})$ (IPP)
$(\rho=100,\ \alpha=\frac{1}{\nu})$

Figure 7: Convergence of Anderson accelerated grad-div stabilized AH method
and IPP iteration for $Re=100$ for varying parameters and $\gamma=100$.

##  6 Conclusions

This paper developed multiple improvements to the AH method for solving the
steady Navier-Stokes equations, and showed that with grad-div stabilization,
SV elements and Anderson acceleration, the AH method can be a very effective
and efficient solver. SV elements and grad-div stabilization allowed us to
connect AH to the well known iterated penalty Picard method, which has good
convergence properties under small data [33]. We also proved that the AH
iteration, under certain conditions on the data and parameters, fits into the
Anderson acceleration analysis framework developed in [30] and thus AA
improves the linear convergence rate of the AH method by the gain of the
underlying AA optimization problem. We also gave results of several numerical
tests that show how each of these improvements is important for good
convergence behavior, and when used together the AH method can be very
effective.

##  7 Acknowledgment

Author PG acknowledges partial support from National Science Foundation grant
DMS 1907823. Authors LR and DV acknowledge partial support from NSF grant DMS
2011490.

## References

  * [1] H. An, X. Jia, and H. Walker.  Anderson acceleration and application to the three-temperature energy equations.  Journal of Computational Physics, 347:1–19, 2017.
  * [2] D. G. Anderson.  Iterative procedures for nonlinear integral equations.  J. Assoc. Comput. Mach., 12(4):547–560, 1965.
  * [3] K. J. Arrow and L. Hurwicz.  Gradient method for concave programming I: Local results.  In L. H. K. J. Arrow and H. Uzawa, editors, Studies in Linear and Nonlinear Programming, page 117–126. Stanford University Press, Stanford, CA, 1958.
  * [4] M. Benzi, G. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numerica, pages 1–137, 2005.
  * [5] M. Benzi and M. Olshanskii.  An augmented Lagrangian-based approach to the Oseen problem.  SIAM J. Sci. Comput., 28:2095–2113, 2006.
  * [6] S. Börm and S. Le Borne.  $\mathcal{H}$-LU factorization in preconditioners for augmented Lagrangian and grad-div stabilized saddle point systems.  Internat. J. Numer. Methods Fluids, 68(1):83–98, 2012.
  * [7] P. Chen, J. Huang, and H. Sheng.  Solving steady incompressible Navier-Stokes equations by the Arrow-Hurwicz method.  Journal of Computational and Applied Mathematics, 311:100–114, 2017\.
  * [8] P. Chen, J. H. P, and H. Sheng.  Some Uzawa methods for steady incompressible Navier-Stokes equations discretized by mixed element methods.  J. Comput. Appl. Math., 273:313–325, 2015.
  * [9] R. Codina.  An iterative penalty method for the finite element solution of the stationary Navier-Stokes equations.  Computer Methods in Applied Mechanics and Engineering, 110:237–262, 1993.
  * [10] C. Evans, S. Pollock, L. Rebholz, and M. Xiao.  A proof that Anderson acceleration improves the convergence rate in linearly converging fixed-point methods (but not in those converging quadratically).  SIAM Journal on Numerical Analysis, 58:788–810, 2020.
  * [11] L. Franca and T. Hughes.  Two classes of mixed finite element methods.  Computer Methods in Applied Mechanics and Engineering, 69(1):89–129, 1988.
  * [12] A. Fu, J. Zhang, and S. Boyd.  Anderson accelerated Douglas-Rachford splitting.  SIAM Journal on Scientific Computing, 42(6):A3560–A3583, 2020.
  * [13] U. Ghia, K. N. Ghia, and C. T. Shin.  High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method.  J. Comput. Phys., 48:387–411, 1982.
  * [14] V. Girault and P.-A. Raviart.  Finite element methods for Navier-Stokes equations, volume 5 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1986.  Theory and algorithms.
  * [15] P. Gresho and R. Lee.  Don’t suppress the wiggles. - they’re telling you something.  Computers and Fluids, 9:223 –253, 1981.
  * [16] M. Gunzburger.  Finite element methods for viscous incompressible flows: A guide to theory, practice, and algorithm (Computer Science and Scientific Computing.  Academic Press Inc.,Boston,MA, 1989.
  * [17] M. Gunzburger.  Iterative penalty methods for the Stokes and Navier-Stokes equations.  Proceedings from Finite Element Analysis in Fluids conference, University of Alabama, Huntsville, pages 1040–1045, 1989.
  * [18] J. Guzman and L. Scott.  The Scott-Vogelius finite elements revisited.  Math. Comp., 88(316):515–529, 2019.
  * [19] T. Heister and G. Rapin.  Efficient augmented Lagrangian-type preconditioning for the Oseen problem using grad-div stabilization.  Int. J. Numer. Meth. Fluids, 71:118–134, 2013.
  * [20] N. Higham and N. Strabic.  Anderson acceleration of the alternating projections method for computing the nearest correlation matrix.  Numerical Algorithms, 72:1021–1042, 2016.
  * [21] V. John and A. Liakos.  Time dependent flow across a step: the slip with friction boundary condition.  International Journal for Numerical Methods in Fluids, 50:713–731, 2006.
  * [22] V. John, A. Linke, C. Merdon, M. Neilan, and L. G. Rebholz.  On the divergence constraint in mixed finite element methods for incompressible flows.  SIAM Review, 59(3):492–544, 2017.
  * [23] C. Kelley.  Numerical methods for nonlinear equations.  Acta Numerica, 27:207–287, 2018.
  * [24] W. Layton.  An Introduction to the Numerical Analysis of Viscous Incompressible Flows.  SIAM, Philadelphia, 2008.
  * [25] J. Loffeld and C. Woodward.  Considerations on the implementation and use of Anderson acceleration on distributed memory and GPU-based parallel computers.  Advances in the Mathematical Sciences, pages 417–436, 2016.
  * [26] P. A. Lott, H. F. Walker, C. S. Woodward, and U. M. Yang.  An accelerated Picard method for nonlinear systems related to variably saturated flow.  Adv. Water Resour., 38:92–101, 2012.
  * [27] H. Morgan and L. Scott.  Towards a unified finite element method for the stokes equations.  SIAM Journal on Scientific Computing, 40(1):A130–A141, 2018.
  * [28] M. A. Olshanskii and A. Reusken.  Grad-Div stabilization for the Stokes equations.  Math. Comp., 73:1699–1718, 2004.
  * [29] Y. Peng, B. Deng, J. Zhang, F. Geng, W. Qin, and L. Liu.  Anderson acceleration for geometry optimization and physics simulation.  ACM Transactions on Graphics, 42:1–14, 2018.
  * [30] S. Pollock and L. Rebholz.  Anderson acceleration for contractive and noncontractive operators.  IMA Journal of Numerical Analysis, 41(4):2841–2872, 2021.
  * [31] S. Pollock, L. Rebholz, and M. Xiao.  Anderson-accelerated convergence of Picard iterations for incompressible Navier-Stokes equations.  SIAM Journal on Numerical Analysis, 57:615– 637, 2019.
  * [32] S. Pollock, L. Rebholz, and M. Xiao.  Acceleration of nonlinear solvers for natural convection problems.  Journal of Numerical Mathematics, 29(4):323–341, 2021.
  * [33] L. Rebholz, D. Vargun, and M. Xiao.  Enabling fast convergence of the iterated penalty Picard iteration with $o(1)$ penalty parameter for incompressible Navier-Stokes via Anderson acceleration.  Computer Methods in Applied Mechanics and Engineering, 387(114178):1–17, 2021.
  * [34] L. Rebholz, A. Viguerie, and M. Xiao.  Efficient nonlinear iteration schemes based on algebraic splitting for the incompressible Navier-Stokes equations.  Math. Comp., 88:1533–1557, 2019.
  * [35] L. Rebholz and M. Xiao.  On reducing the splitting error in Yosida methods for the Navier-Stokes equations with grad-div stabilization.  Computer Methods in Applied Mechanics and Engineering, 294:259–277, 2015.
  * [36] P. Stasiak and M. Matsen.  Efficiency of pseudo-spectral algorithms with anderson mixing for the SCFT of periodic block-copolymer phases.  Eur. Phys. J. E, 34:110:1–9, 2011.
  * [37] R. Temam.  Navier-Stokes equations. Theory and numerical analysis.  North-Holland Publishing Co., Amsterdam, 1977.  Studies in Mathematics and its Applications, Vol. 2.
  * [38] A. Toth, C. Kelley, S. Slattery, S. Hamilton, K. Clarno, and R. Pawlowski.  Analysis of Anderson acceleration on a simplified neutronics/thermal hydraulics system.  Proceedings of the ANS MC2015 Joint International Conference on Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method, ANS MC2015 CD:1–12, 2015.
  * [39] H. Uzawa.  Iterative methods for concave programming.  In L. H. K. J. Arrow and H. Uzawa, editors, Studies in Linear and Nonlinear Programming, page 154–165. Stanford University Press, Stanford, CA, 1958.
  * [40] H. F. Walker and P. Ni.  Anderson acceleration for fixed-point iterations.  SIAM J. Numer. Anal., 49(4):1715–1735, 2011.
  * [41] D. Wicht, M. Schneider, and T. Bohlke.  Anderson‐accelerated polarization schemes for fast Fourier transform‐based computational homogenization.  International Journal for Numerical Methods in Engineering, to appear, 2021.