# Efficient Solvers for Large-Scale Saddle Point Systems Arising in Feedback Stabilization of Multi-field Flow Problems

Peter Benner, Jens Saak, Martin Stoll, Heiko K. Weichelt

## ► To cite this version:

Peter Benner, Jens Saak, Martin Stoll, Heiko K. Weichelt. Efficient Solvers for Large-Scale Saddle Point Systems Arising in Feedback Stabilization of Multi-field Flow Problems. 26th Conference on System Modeling and Optimization (CSMO), Sep 2013, Klagenfurt, Austria. pp.11-20, 10.1007/978-3-662-45504-3_2. hal-01286209

HAL Id: hal-01286209

https://inria.hal.science/hal-01286209

Submitted on 10 Mar 2016

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Efficient solvers for large-scale saddle point
systems arising in feedback stabilization of
multi-field flow problems

Peter Benner, Jens Saak, Martin Stoll, and Heiko K. Weichelt

Max Planck Institute for Dynamics of Complex Technical Systems

Sandtorstr. 1, 39106 Magdeburg, Germany

`{benner,saak,stoll,weichelt}@mpi-magdeburg.mpg.de`

http://www.mpi-magdeburg.mpg.de

**Abstract.** This article introduces a block preconditioner to solve large-scale block structured saddle point systems using a Krylov-based method. Such saddle point systems arise, e.g., in the Riccati-based feedback stabilization approach for multi-field flow problems as discussed in [2]. Combining well known approximation methods like a least-squares commutator approach for the Navier-Stokes Schur complement, an algebraic multigrid method, and a Chebyshev-Semi-Iteration, an efficient preconditioner is derived and tested for different parameter sets by using a simplified reactor model that describes the spread concentration of a reactive species forced by an incompressible velocity field.

**Keywords:** Coupled flow control, large-scale saddle point systems, pre-conditioned GMRES, least-squares commutator approach, algebraic multi-grid, Chebyshev-Semi-Iteration

# 1 Introduction

In this paper we investigate the solution of large-scale saddle point systems arising in control problems for coupled partial differential equations (PDEs). The starting points are recent publications concerning the boundary feedback stabilization of non-coupled flows like the linear Stokes flow in [3] and the non-linear Navier-Stokes flow in [1]. The analytic approach for this feedback stabilization is given by Raymond in, e.g., [13].

Using the projection idea proposed by Heinkenschloss et al. [8], Benner et al. [1,3] show that the solution of certain saddle point systems is the key ingredient to ensure that the numerical solution lies on the correct solution manifold, i.e., the space of discretely divergence-free velocity fields, without performing an explicit projection.

Applying these ideas to a coupled flow problem, namely the Navier-Stokes equations combined with a diffusion-convection equation, leads to saddle point systems with a more complicated block structure [2]. Solving these systems efficiently requires the use of appropriate preconditioners. This paper investigates
an efficient iterative solution strategy via the use of preconditioned Krylov sub-space methods based on the framework derived in [3]. Here we consider the full feedback system for the coupled multi-field flow problem, while in [3], only the linear Stokes case was treated without coupling to another field equation. Moreover, this paper complements [2] in the sense that there, we have focused on presenting results on the convergence of the Newton-ADI method for solving the algebraic Riccati equation determining the stabilizing feedback control for the coupled system, where the saddle point problems in the innermost step of the Newton-ADI iteration were solved by sparse direct methods, while here, we study preconditioned iterative solvers for this step.

The remainder of this paper is organized as follows. Section 2 briefly recalls the feedback stabilization approach for multi-field flow problems from [2] that leads to large-scale saddle point systems. Afterwards, we discuss properties of these saddle point systems to derive a suitable preconditioner in Section 3. Section 4 shows numerical results before we conclude the paper and give a short outlook to further investigations in Section 5.

## 2 Derivation of saddle point systems

The derivation of the block structured saddle point systems in [2] starts with the linearized coupled flow problem defined for $t \in [0, \infty)$ and $x \in \Omega \subset \mathbb{R}^2$. The linearized Navier-Stokes equations that describe, up to first order, the difference between actual and desired velocity and pressure are given as

$$ \frac{\partial z}{\partial t} - \frac{1}{\text{Re}} \Delta z + (w \cdot \nabla)z + (z \cdot \nabla)w + \nabla p = f_l, \quad \text{on } [0, \infty) \times \Omega. \quad (1) \\ \text{div } z = 0, $$

They are then coupled via the velocity field $z(t, x)$ with the linearized diffusion-convection equation

$$ \frac{\partial c_z}{\partial t} - \frac{1}{\text{ReSc}} \Delta c_z + (w \cdot \nabla)c_z + (z \cdot \nabla)c_w = 0, \quad \text{on } [0, \infty) \times \Omega \quad (2) $$

that describes the concentration of a reactive species denoted by $c_z(t, x)$. The stationary linearization points $w(x)$ for the velocity and $c_w(x)$ for the concentration are assumed to be given. The equations are scaled with the Reynolds number $Re$ and the Schmidt number $Sc$. Using the mixed Taylor-Hood finite elements [9] for the velocity and pressure in Eqs. (1) as well as linear ansatz functions for the concentration in Eq. (2), we end up with a system of discrete differential-algebraic equations (DAE) that can be written as the control system:
$$ \begin{bmatrix} M_z & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & M_c \end{bmatrix} \frac{d}{dt} \begin{bmatrix} \mathbf{z} \\ \mathbf{p} \\ \mathbf{c} \end{bmatrix} = \begin{bmatrix} A_z & G & 0 \\ G^T & 0 & 0 \\ -R & 0 & A_c \end{bmatrix} \begin{bmatrix} \mathbf{z} \\ \mathbf{p} \\ \mathbf{c} \end{bmatrix} + \begin{bmatrix} B_z \\ 0 \\ 0 \end{bmatrix} \mathbf{u}, \quad (3a) $$

$$ \mathbf{y} = \begin{bmatrix} 0 & 0 & C_c \end{bmatrix} \begin{bmatrix} \mathbf{z} \\ \mathbf{p} \\ \mathbf{c} \end{bmatrix} \quad (3b) $$

with the first block row for velocity (of dimension $n_z$), the second row for pressure (of dimension $n_p$), and the third row for concentration (of dimension $n_c$) [2].
The matrix pencil

$$ \underbrace{\begin{bmatrix} A_z & G & 0 \\ G^T & 0 & 0 \\ -R & 0 & A_c \end{bmatrix}}_{\mathbf{A}}, \underbrace{\begin{bmatrix} M_z & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & M_c \end{bmatrix}}_{\mathbf{M}} $$

is of dimension $n \times n$ with $n = n_z + n_c + n_p$ and has $2n_p$ infinite eigenvalues [5].

In [2], a linear-quadratic regulator (LQR) approach is applied to system (3) for determining the stabilizing control function $\mathbf{u}$. The solution of this LQR problem is a linear feedback control $\mathbf{u}(t) = \mathbf{K}(\mathbf{z}(t), \mathbf{p}(t), \mathbf{c}(t))$, determined via the solution of an algebraic Riccati equation (ARE) defined on the subspace of discretely divergence-free vector fields. The resulting ARE is then solved using a Newton-ADI algorithm. This method yields a threefold nested iteration. In the innermost loop, saddle point systems of the form

$$ \underbrace{\begin{bmatrix} A_z^T + q_i M_z & G & -R^T \\ G^T & 0 & 0 \\ 0 & 0 & A_c^T + q_i M_c \end{bmatrix}}_{=A^T+q_i M=:\mathbf{F}_i} \underbrace{\begin{bmatrix} \Lambda_z \\ \Lambda_p \\ \Lambda_c \end{bmatrix}}_{\mathbf{\Lambda}} = \underbrace{\begin{bmatrix} \tilde{Y}_z \\ 0 \\ \tilde{Y}_c \end{bmatrix}}_{\mathbf{Y}}. \quad (4) $$

have to be solved for certain ADI shifts $q_i \in \mathbb{C}^-$ and a block right hand side $\mathbf{Y}$. The whole nested iteration is given in [2, Algorithm 1] and is omitted here due to space constraints.

### 3 Preconditioned iterative solvers for block structured saddle point systems

The use of direct solvers in (4) is only suitable for moderate problem sizes and two-dimensional problems. Although iterative methods can handle much larger systems, their performance will deteriorate if the mesh-size decreases. To avoid this, a suitable preconditioner $\mathbf{P}_i \in \mathbb{C}^{n \times n}$ is introduced such that

$$ \mathbf{P}_i^{-1} \mathbf{F}_i \mathbf{\Lambda} = \mathbf{P}_i^{-1} \mathbf{Y} $$

is solved instead of (4) (see [7,16]). Before we derive a suitable preconditioner $\mathbf{P}_i$ we need to describe the properties of the saddle point system and their influence on the chosen preconditioner.
## 3.1 Properties

The matrices $M_z, M_c$ are symmetric and positive definite, $G, R$ are of full rank, and the ADI shift $q_i \in \mathbb{C}^-$ is contained in the convex hull of the finite spectrum of $(A, M)$. The shifted system matrix $F_i$ is indefinite $\forall q_i \in \mathbb{C}^-$. Due to the different $q_i$, the matrix $F_i$ changes in each ADI step and, therefore, the preconditioner has to be adapted in each ADI step as well. Nevertheless, for the remainder of this section we assume a fixed ADI shift $q_i = q$ to omit the index $i$ if it is obvious.

## 3.2 Derivation of block preconditioner

Adapting the ideas from [3, Section 3.2] we consider

$$ \mathbf{F} = \begin{bmatrix} F_z & G & -R^T \\ G^T & 0 & 0 \\ 0 & 0 & F_c \end{bmatrix} = \begin{bmatrix} \mathbf{F}_{NSE} & -\tilde{R}^T \\ 0 & F_c \end{bmatrix} \quad \text{with} \quad \begin{array}{l} F_z := A_z^T + qM_z, \\ F_c := A_c^T + qM_c, \\ \tilde{R} := [R \ 0], \end{array} \qquad (5) $$

and $\mathbf{F}_{NSE}$ as the saddle point system for the non-coupled Navier-Stokes flow as it is used in [1]. Using the preconditioner $\mathbf{P}_{NSE}$ from [3], we define a block preconditioner for the use with GMRES [17] to solve with the block structured saddle point system (5) as follows:

$$ \mathbf{P} = \begin{bmatrix} \mathbf{P}_{NSE} & -\tilde{R}^T \\ 0 & P_c \end{bmatrix} = \begin{bmatrix} P_z & 0 & -R^T \\ G^T & -P_{SC} & 0 \\ 0 & 0 & P_c \end{bmatrix} \\ \Rightarrow \mathbf{P}^{-1} = \begin{bmatrix} P_z^{-1} & 0 & P_z^{-1} R^T P_c^{-1} \\ P_{SC}^{-1} G^T P_z^{-1} & -P_{SC}^{-1} & P_{SC}^{-1} G^T P_z^{-1} R^T P_c^{-1} \\ 0 & 0 & P_c^{-1} \end{bmatrix}. $$

In contrast to the preconditioner derived in [3], we cannot achieve a block lower triangular matrix due to the coupling matrix $R$. Applying $\mathbf{P}^{-1}$ to $\mathbf{F}$ yields

$$ \mathbf{P}^{-1}\mathbf{F} = \begin{bmatrix} P_z^{-1}F_z & P_z^{-1}G & -P_z^{-1}R^T + P_z^{-1}R^T P_c^{-1}F_c \\ P_{SC}^{-1}G^T P_z^{-1}F_z - P_{SC}^{-1}G^T & P_{SC}^{-1}G^T P_z^{-1}G & -P_{SC}^{-1}G^T P_z^{-1}R^T + P_{SC}^{-1}G^T P_z^{-1}R^T P_c^{-1}F_c \\ 0 & 0 & P_c^{-1}F_c \end{bmatrix} \quad (6) $$

If one assumes $P_z = F_z, P_c = F_c$, and $P_{SC} = G^T F_z^{-1} G$ as ideal approximations in (6), this leads to

$$ \begin{bmatrix} I_z & F_z^{-1}G & -F_z^{-1}R^T + F_z^{-1}R^T \\ P_{SC}^{-1}G^T - P_{SC}^{-1}G^T & P_{SC}^{-1}G^T F_z^{-1}G & -P_{SC}^{-1}G^T F_z^{-1}R^T + P_{SC}^{-1}G^T F_z^{-1}R^T \\ 0 & 0 & I_c \end{bmatrix} = \begin{bmatrix} I_z * 0 \\ 0 & I_p 0 \\ 0 & 0 & I_c \end{bmatrix} $$
and our iterative method would converge within one step. The goal is to find good approximations for $P_z$, $P_c$, and $P_{SC}$ that can be evaluated fast and still cluster the eigenvalues in a suitable way such that our iterative solver shows fast convergence [7]. Instead of calculating the inverse $\mathbf{P}^{-1}$ to apply the preconditioner $\mathbf{P}$, we consider the solution of a linear system

$$ \begin{bmatrix} P_z & 0 & -R^T \\ G^T - P_{SC} & 0 & 0 \\ 0 & 0 & P_c \end{bmatrix} \begin{bmatrix} x_z \\ x_p \\ x_c \end{bmatrix} = \begin{bmatrix} b_z \\ b_p \\ b_c \end{bmatrix} \qquad (7) $$

that can be solved in three steps:

$$ \text{Step I:} \quad x_c = P_c^{-1} b_c, \qquad (8a) $$

$$ \text{Step II:} \quad x_z = P_z^{-1}(R^T x_c + b_z), \qquad (8b) $$

$$ \text{Step III:} \quad x_p = P_{SC}^{-1}(G^T x_z - b_p). \qquad (8c) $$

In conclusion, the coupling matrix $R$ only leads to a matrix-vector multiplication. In steps I and II, one needs to solve with the shifted velocity and concentration system matrices as defined in (5). For both steps, an algebraic multigrid (AMG) method can be used as it is described below. But first, we discuss the more challenging step III that is handled as follows.

## 3.3 Approximation methods

*Schur Complement Approximation.* $P_{SC}$ is an approximation of the Navier-Stokes Schur complement $SC := G^T F_z^{-1} G \in \mathbb{R}^{n_p \times n_p}$. Unfortunately, the matrix $SC$ would be a dense matrix that includes the inverse of $F_z$. To avoid the use of this matrix, we follow the approach in [3,18] and use a slightly modified variant of the least squares commutator approach as it is described in [7, Section 8.2]. Namely, we consider the shifted Oseen operator in the velocity space

$$ \mathcal{F}_z = -\frac{1}{\text{Re}} \nabla^2 + \boldsymbol{w} \cdot \nabla + q\mathcal{I}. $$

Note that it is common practice to omit the reaction term $(\boldsymbol{z} \cdot \nabla)\boldsymbol{w}$ that appears in the linearized Navier-Stokes equations to derive preconditioners [7, Section 8]. Similar to [7, Section 8.2] and [6], we suppose that there exists an analogous operator on the pressure space defined as

$$ \mathcal{F}_p = (-\frac{1}{\text{Re}} \nabla^2 + \boldsymbol{w} \cdot \nabla + q\mathcal{I})_p. $$

The least squares commutator of the shifted Oseen operator with the gradient operator is defined as

$$ \mathcal{E} = (\mathcal{F})\nabla - \nabla(\mathcal{F}_p) $$
and is supposed to become small in some sense [7]. Using the discrete versions of the operators, we end up with

$$E = (M_{\mathbf{z}}^{-1} F_{\mathbf{z}}) M_{\mathbf{z}}^{-1} G - M_{\mathbf{z}}^{-1} G (M_{\mathbf{p}}^{-1} F_{\mathbf{p}})$$

with $M_{\mathbf{p}}$ the mass matrix and $F_{\mathbf{p}} = A_{\mathbf{p}}^{T} + qM_{\mathbf{p}}$ the shifted system matrix, both defined on the pressure space. Premultiplying this by $G^{T}F_{\mathbf{z}}^{-1}M_{\mathbf{z}}^{-1}$ and postmultiplying by $F_{\mathbf{p}}^{-1}M_{\mathbf{p}}$ yields

$$G^{T} M_{\mathbf{z}}^{-1} G F_{\mathbf{p}}^{-1} M_{\mathbf{p}} \approx G^{T} F_{\mathbf{z}}^{-1} G = SC,$$

[3]. The large and dense matrix $G^{T}M_{\mathbf{z}}^{-1}G$ cannot be used explicitly, but it is shown in [7, Section 5.5.1] that this matrix is spectrally equivalent to the Laplacian $S_{\mathbf{p}}$ defined on the pressure space for an inf-sup stable discretization and an inflow-outflow problem [7, Section 8.2], as it is considered in this paper. Finally, we obtain

$$P_{SC} \approx S_{\mathbf{p}} F_{\mathbf{p}}^{-1} M_{\mathbf{p}} \Rightarrow P_{SC}^{-1} \approx M_{\mathbf{p}}^{-1} F_{\mathbf{p}} S_{\mathbf{p}}^{-1}.$$

In [4] the authors use a similar approach for the Navier-Stokes equations. In summary, the application of $P_{SC}^{-1}$ requires to solve with $S_{\mathbf{p}}$ (step IIIa), multiply with $F_{\mathbf{p}}$ (step IIIb), and solve with $M_{\mathbf{p}}$ (step IIIc). The step IIIa can be solved with an AMG method, similar to the steps I and II.

*Algebraic Multigrid.* As it is described above, the steps I (8a), II (8b) and IIIa are solved using an AMG method [14]. Due to the possibly complex ADI shifts $q$ in (8a) and (8b), we use the AGMG package developed by the group of Y. Notay [10,11,12]. In all three cases we use the MATLAB®-based implementation to solve systems of the form

$$Fx = b$$

with a sparse matrix $F \in \{F_{\mathbf{z}}, F_{\mathbf{c}}, S_{\mathbf{p}}\}$. Details about the used parameters for the function `agmg` are discussed in Subsection 4.2. For more details about the internally used methods and the implemented syntax we refer the reader to [11]. Although the AGMG method can handle complex arithmetic, it needs significantly more steps to converge to the desired tolerance. Additionally, we note that `agmg` is a non-linear function, such that one should use a flexible iterative method, e.g., FGMRES [15]. However, our numerical experiments do not show any drawbacks using a standard GMRES implementation.

*Chebyshev-Semi-Iteration.* Although the solution of step IIIc with the symmetric positive definite mass matrix $M_{\mathbf{p}}$ is relatively cheap, this can still be accelerated by using the Chebyshev-Semi-Iteration as it is described, e.g., in [18]. Numerical tests showed that one needs only 4–6 steps to obtain a suitable result for the preconditioner, which results in a speedup that is shown in Subsection 4.2.

The next section depicts selected results to show the performance of the preconditioned iterative method.
<figure><img src="image_2.png" /><figcaption>Fig. 1: Initial triangulation of the reactor model with coordinates and boundary conditions [2].</figcaption></figure>

<table><thead><tr><th>Set</th><th>Re</th><th>Sc</th></tr></thead><tbody><tr><td>I</td><td>1</td><td>1</td></tr><tr><td>II</td><td>1</td><td>10</td></tr><tr><td>III</td><td>10</td><td>1</td></tr><tr><td>IV</td><td>1</td><td>100</td></tr><tr><td>V</td><td>10</td><td>10</td></tr></tbody></table>

(a) Different parameter settings.

<table><thead><tr><th>Variable</th><th>Dimension</th></tr></thead><tbody><tr><td>$n_z$</td><td>9 092</td></tr><tr><td>$n_p$</td><td>1 276</td></tr><tr><td>$n_c$</td><td>1 187</td></tr><tr><td>$n$</td><td>11 555</td></tr></tbody></table>

(b) Different dimensions of FE space.

<caption>Table 1: Test parameter settings.</caption>

# 4 Numerical examples

To test the efficiency of the preconditioned iterative method, the same data and configurations as in [2] are used. After refining the initial triangulation of the reactor model in Figure 1, we end up with the variable dimensions as depicted in Table 1b. Furthermore, we define five parameter sets for different combinations of Reynolds and Schmidt numbers as shown in Table 1a. We use the MATLAB implementation of GMRES [17] to solve the saddle point systems (4) for selected ADI shifts $q_i$ that appear during the Newton-ADI iteration. Each $q_i$ is used for three ADI steps with four right hand sides every time. Thereby, the number of GMRES steps and the CPU times are measured and arithmetically averaged. The preconditioner $P$ is evaluated as a MATLAB function handle that solves the linear system (7) using the steps (8). The GMRES tolerance is set to $10^{-10}$ to ensure the same convergence of the ADI iteration that a direct solve would imply [3]. Although a few complex ADI-shifts $q_i$ appear for each parameter set during the Newton-ADI process, the pictures only show the real parts of $q_i$.

All computations were executed in MATLAB R2012a on a 64-bit server with 2×Intel® Xeon® X5650 @2.67GHz, 12 Cores (6 Cores per CPU) and 48 GB main memory available.

## 4.1 Influence of ADI shifts and Reynolds and Schmidt numbers

The influence of the variation of the Reynolds and Schmidt numbers as given in Table 1 is depicted in Figure 2. To obtain the best approximations for the
<figure><img src="image_3.png" /><figcaption>Fig. 2: Average number of GMRES steps for a representative selection of ADI shifts from the Newton-ADI iteration for the configuration sets in Table 1a.</figcaption></figure>

preconditioning steps (8a) – (8c), a direct solver is used to solve with $F_z, F_c$, and $S_p$. It can be observed that for ADI shifts $-10^5 < \text{Re}(q_i) < -10^1$, between 20 – 25 GMRES steps are needed. As soon as the absolute value of $q_i$ gets smaller then 10 the number of steps increases. This is a natural behavior, because the influence of the mass matrices $M_z$ and $M_p$ vanishes. Nevertheless, GMRES converges within at most 40 – 80 steps for all parameter configurations. An empirical test to set: $q_i = -10 \quad \forall |q_i| < 10$, during the Newton-ADI process showed similar ADI convergence behavior as for the original shift selection, without the drawback of higher GMRES cost for certain shifts. In summary, the derived preconditioner is suitable concerning different Reynolds and Schmidt numbers, as well as different ADI shifts.

## 4.2 Approximations using AMG and Chebyshev-semi-iteration

As described in Subsection 3.3, the different preconditioning steps should be solved by an easy to evaluate approximation that is accurate enough to ensure the convergence of GMRES, but avoids the use of sparse factorizations of large-scale sparse matrices. We exchanged the direct solver by its approximation step by step and depict the results in Figure 3. At first, we use the MATLAB based function `agmg` [11] to solve with $F_z$ and $F_c$ in (8b) and (8a) with an accuracy of $10^{-10}$. Depending on the used ADI shift, the function `agmg` needed 1 – 30 steps. Thus, the times to solve the whole saddle point system with the same number of GMRES steps increased a little bit compared to the direct solver. At second, we approximately solved with $S_p$ in step IIIa using `agmg` as a preconditioner. This was sufficient enough to achieve the GMRES accuracy and, furthermore, decreased the time. Finally, we applied a Chebyshev-Semi-Iteration to approximately solve with $M_p$ in step IIIc. The obtained speedup finally decreased the times below the time used by the direct solver in each step without the loss of
<figure><img src="image_4.png" /><figcaption>Fig. 3: Average time to solve Eq. (4) with GMRES for a representative selection of real ADI shifts from the Newton-ADI iteration for different approximations of the preconditioning steps (8).</figcaption></figure>

any accuracy in GMRES. Due to the above addressed problems with complex ADI shifts in agmg, we restrict our comparison in Figure 3 to a selection of real ADI shifts. The selection has been performed such that the span of all ADI shifts appearing in the entire Newton-ADI process is covered. Where those shifts clustered we only chose one representative per cluster.

At the end of this section it should be noted that the suggested preconditioned GMRES method for the considered class of saddle point problems would show its full strength in comparison to a direct solver when using finer discretizations, leading to larger dimensions, and in particular when moving to 3D problems. This will be addressed in future work.

# 5 Conclusions and outlook

We have recalled the formation of block structured saddle point systems as they arise within the Riccati-based feedback stabilization approach for coupled flow problems that avoids any explicit projection [2]. We were able to extend the results from [3], developed for uncoupled Stokes flow, to the coupled flow described by incompressible Navier-Stokes and a diffusion-convection equation. For that reason, the least-squares commutator approach in [7] has been modified to approximate the shifted Navier-Stokes Schur complement. Exploiting the block structure of the arising preconditioner guarantees a fast evaluation within GMRES. Each of the blocks can either be approximated by an AMG method or a Chebyshev-Semi-Iteration. Several numerical experiments showed that the derived preconditioning method is able to solve the arising saddle point systems efficiently independent of the different parameter settings. Only the use of complex ADI shifts during the Newton-ADI process is not yet optimally covered by this approach and will be investigated in the future.
References

1. Bänsch, E., Benner, P., Saak, J., Weichelt, H.K.: Riccati-based boundary feedback stabilization of incompressible Navier-Stokes flow. Preprint SPP1253-154, DFG-SPP1253 (2013)

2. Bänsch, E., Benner, P., Saak, J., Weichelt, H.K.: Optimal control-based feedback stabilization of multi-field flow problems. In: Leugering, G., Benner, P., Engell, S., Griewank, A., Harbrecht, H., Hinze, M., Rannacher, R., Ulbrich, S. (eds.) Trends in PDE Constrained Optimization. International Series of Numerical Mathematics, Birkhäuser (2014, to appear)

3. Benner, P., Saak, J., Stoll, M., Weichelt, H.K.: Efficient solution of large-scale saddle point systems arising in Riccati-based boundary feedback stabilization of incompressible Stokes flow. SIAM J. Sci. Comput. 35(5), S150–S170 (2013)

4. Benzi, M., Olshanskii, M.A., Wang, Z.: Modified augmented Lagrangian preconditioners for the incompressible Navier-Stokes equations. Int. J. Numer. Methods Fluids 66(4), 486–508 (2011)

5. Cliffe, K.A., Garratt, T.J., Spence, A.: Eigenvalues of block matrices arising from problems in fluid mechanics. SIAM J. Matrix Anal. Appl. 15(4), 1310–1318 (1994)

6. Elman, H., Howle, V., Shadid, J., Shuttleworth, R., Tuminaro, R.: Block preconditioners based on approximate commutators. SIAM J. Sci. Comput. 27(5), 1651–1668 (2006)

7. Elman, H., Silvester, D., Wathen, A.: Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics. Oxford University Press, Oxford (2005)

8. Heinkenschloss, M., Sorensen, D.C., Sun, K.: Balanced truncation model reduction for a class of descriptor systems with application to the Oseen equations. SIAM J. Sci. Comput. 30(2), 1038–1063 (2008)

9. Hood, P., Taylor, C.: Navier-Stokes equations using mixed interpolation. In: Oden, J.T., Gallagher, R.H., Taylor, C., Zienkiewicz, O.C. (eds.) Finite Element Methods in Flow Problems. pp. 121–132. University of Alabama in Huntsville Press (1974)

10. Napov, A., Notay, Y.: An algebraic multigrid method with guaranteed convergence rate. SIAM J. Sci. Comput. 34, A1079–A1109 (2012)

11. Notay, Y.: An aggregation-based algebraic multigrid method. ETNA, Electron. Trans. Numer. Anal. 37, 123–146 (2010)

12. Notay, Y.: Aggregation-based algebraic multigrid for convection-diffusion equations. SIAM J. Sci. Comput. 34, A2288–A2316 (2012)

13. Raymond, J.P.: Feedback boundary stabilization of the two-dimensional Navier-Stokes equations. SIAM J. Control Optim. 45(3), 790–828 (2006)

14. Ruge, J., Stüben, K.: Algebraic multigrid (amg). In: S. McCormich, E. (ed.) Multi-grid Methods. Frontiers in Applied Mathematics, vol. 5, pp. 73–130. SIAM (1987)

15. Saad, Y.: A flexible inner-outer preconditioned gmres algorithm. SIAM J. Sci. Comput. 14(2), 461–469 (1993)

16. Saad, Y.: Iterative methods for sparse linear systems. SIAM, Philadelphia, PA, USA (2003)

17. Saad, Y., Schultz, M.H.: GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. and Stat. Comput. 7(3), 856–869 (1986)

18. Stoll, M., Wathen, A.: All-at-once solution of time-dependent Stokes control. J. Comp. Phys. 232, 498–515 (2013)