A short course on numerical simulation of viscous flow:Discretization, optimization and stability analysis

This article contains part of the material of four introductory lectures given 
at the 12th school ``Mathematical Theory in Fluid Mechanics'', Spring 2011, 
at Kacov, Czech Republic, on ``Numerical simulation of viscous flow: 
discretization, optimization and stability analysis''. 
In the first lecture on ``Numerical computation of incompressible viscous flow'', 
we discuss the Galerkin finite element method for the discretization of the 
Navier-Stokes equations for modeling laminar flow. Particular emphasis 
is put on the aspects pressure 
stabilization and truncation to bounded domains. 
In the second lecture on ``Goal-oriented adaptivity'', we introduce the concept 
underlying the Dual Weighted Residual (DWR) method for goal-oriented residual-based 
adaptivity in solving the Navier-Stokes equations. This approach is 
presented for stationary as well as nonstationary situations. 
In the third lecture on ``Optimal flow control'', we discuss the use of the 
DWR method for adaptive discretization in flow control and model calibration. 
Finally, in the fourth lecture on ``Numerical stability analysis'', we consider 
the numerical stability analysis of stationary flows employing the concepts 
of linearized stability and pseudospectrum.

1. Numerical computation of incompressible viscous flow. In this first lecture, we discuss the Galerkin finite element method for the discretization of the Navier-Stokes equations for modeling laminar flow, i.e., for the case that all scales of the flow can be properly resolved. Particular emphasis is put on the aspects pressure stabilization and truncation to bounded domains. For further information, we refer to the relevant text book literature Pironneau [69], Girault & Raviart [41], Quartapelle [71], Fertzinger & Peric [35], Gresho & Sani [44], and to the more recent survey articles Glowinski [42], Rannacher [73,74,75].
A state equation models the connection of pressure and density. In the following, the fluid is assumed to be of Stokes-type, and incompressible, and the density to be homogeneous, ρ ≡ ρ 0 . In this case equation (1) reduces to the constraint ∇ · v = 0 and only the viscosity coefficient µ appears in the momentum equation (2). Then, in the isothermal case the Navier-Stokes system can be written in short as with the kinematic viscosity ν = µ/ρ 0 > 0. This system is supplemented by initial and boundary conditions v| t=0 = v 0 , v| Γ rigid = 0, v| Γin = v in , µ∂ n v + pn| Γout = 0, where Γ rigid , Γ in , Γ out are the rigid part, the inflow part and the outflow part, respectively, of the flow domain's boundary. The role of the natural outflow boundary condition on Γ out will be discussed in more detail below. In this formulation the flow domain may be two-or three-dimensional. This model is made dimensionless through a scaling transformation with the Reynolds number Re = U L/ν as the characteristic parameter, where U is the reference velocity and L the characteristic length, e.g., U ≈ max |v in | and L ≈ diam(Ω). The Galerkin Finite Element Method (FEM) is based on a variational formulation of the problem and determines "discrete" approximations of the primitive variables velocity and pressure in certain finite dimensional trial spaces consisting of piecewise polynomial functions. This discretization inherits most of the structure of the continuous problem, which results in high computational flexibility and provides the basis of a solid mathematical analysis.
It is well known that in two space dimensions the pure Dirichlet problem (8), (9), with Γ out = ∅, possesses a unique solution on any time interval [0, T ], which is also a classical solution if the data of the problem are regular enough. However, for small viscosity, i.e., large Reynolds number, this solution may be unstable.
In three dimensions, the existence of a unique solution is known only for sufficiently small data, e.g., ∇v 0 ≈ ν, or on sufficiently short intervals of time, 0 ≤ t ≤ T , with T ≈ ν. The proof or disproof of the general result is considered as one of the most challenging theoretical problems in Mathematical Fluid Mechanics.

Remark 1. Sometimes the (vector) Burgers equation
is considered as a simplified model of the Navier-Stokes equation since it involves a nonlinearity of the same form. However, due to the missing incompressibility constraint the Burgers equation has significantly different properties than the Navier-Stokes equations. For example, in contrast to the latter the Burgers equation admits a maximum principle, which allows for the proof of global solvability, and also its spectral properties are very different from those of the Navier-Stokes equations. The latter aspect will be discussed in more detail in the context of hydrodynamic stability, below.
1.2. Regularity of solutions. Formally, even for large Reynolds number, the Navier-Stokes equations are of elliptic type in the stationary and of parabolic type in the nonstationary case. This means that qualitatively the solution of the problem should be as regular as the data and the boundary of the flow domain permit. However, numerical approximation depends on quantitative regularity properties of the solution. For example, the presence of a strong boundary layer or part of the boundary with large curvature may strongly effect the accuracy on realistic meshes, though the solution is qualitatively smooth. The regularity estimates for the solutions of the Navier-Stokes equations provided in the mathematical literature are mostly of such a qualitative character and have only limited value in predicting the performance of a concrete numerical computation. Therefore, we do not list these results but rather concentrate on some aspects of regularity which are peculiar to the incompressible Navier-Stokes equations. For the general existence and regularity theory for the Navier-Stokes equations, we refer to the standard literature, e.g., Ladyshenskaya [60], Heywood [47], Temam [88], Feistauer [34], Galdi [38], Friedlander & Serre [37].
1.2.1. Solution behavior as t → 0. For the purposes of numerical analysis one needs regularity of the solution uniformly down to t = 0, which turns out to be a delicate requirement. To illustrate this, let us assume that the solution {v, p} is uniformly smooth as t → 0 . Then, applying the divergence operator to the Navier-Stokes 1.2.2. Solution behavior as t → ∞. The constants in the a priori error estimates available in the literature for space and time discretization of the Navier-Stokes equation generally depend exponentially on the Reynolds number Re ≈ ν −1 and on the final time T . In general, this exponential dependence is unavoidable, due to the use of Gronwall's inequality in the derivation of these estimates, and may also be natural in the case of unstable solutions. Since even in the range of laminar flows e 20 ≈ 10 8 , e 100 ≈ 10 43 , the practical meaning of a priori error estimates containing such constants is questionable (Johnson et al. [58]). If the solution to be computed were known to possess certain stability properties, for instance to be "exponentially stable", then the exponential dependence on T can be suppressed, so justifying long-time flow simulation (Heywood & Rannacher [49]). However, such stability assumptions are usually not verifiable by a priori analysis but may be justified in many situations by experimental evidence. Consequently, reliable flow simulation requires computable a posteriori error bounds in terms of the approximate solution. This will be the subject of one of the next sections.
1.3. Artificial boundary conditions: truncation to bounded domains. The variational formulation (8), (9) does not contain an explicit reference to any "outflow boundary condition". Suppose that the solution v ∈ v in + H, p ∈ L is sufficiently smooth. Then, integration by parts and varying ϕ ∈ H in the terms ν(∇v, ∇ϕ) − (p, ∇ · ϕ) = (ν∂ n v − pn, ϕ) Γout + (−ν∆v + ∇p, ϕ) yields the already mentioned "natural" outflow boundary condition This condition has proven to be well suited in modeling almost parallel flows (Heywood et al. [53]). It occurs in the variational formulation of the problem if one does not prescribe any boundary condition for the velocity at the outlet suggesting the name "do nothing" or "free" boundary condition. However, this approach has to be used with care. For example, if the flow region contains more than one outlet, undesirable effects may occur, since the "do nothing" condition contains an additional hidden condition on the pressure. In fact, integrating (12) over any component S of the outflow boundary and using the incompressibility constraint ∇ · v = 0 and v |Γ rigid = 0 yields i.e., the mean pressure over S must be zero. To illustrate this, let us consider low Reynolds number flow through a junction in a system of pipes, again prescribing a Poiseuille inflow upstream (see Figure 1). Obviously, making one leg of the pipe longer significantly changes the flow pattern. By the property (13) the pressure gradient is greater in the shorter of the two outflow sections, which explains why there is a greater throughflow. . This change has no effect in the case of pure Dirichlet boundary conditions as then the two forms coincide. But in using the "do nothing" approach this modification leads to the outflow boundary condition which results in outward bending streamlines not representing Poiseuille flow as assumed in an infinitely extended pipe (Heywood et al. [53]).
Remark 3. Despite its remarkable success of the "do nothing" outflow boundary condition in modeling (almost) parallel flow, there is a theoretical problem with existence and uniqueness (Heywood et al. [53], Rannacher [75]). Existence of solutions (in 2d as well as in 3d) can be shown even in the stationary case only for very small Reynolds numbers. The problem is the lacking a priori bound for the Galerkin approximations in the Dirchlet norm due to the indefiniteness of the nonlinear term: This dilemma can also be described in terms of uniqueness: Is v ≡ 0 the only solution to the following homogeneous problem?
Extensive numerical tests did not lead to an affirmative answer to this question. Therefore, the use of the "do nothing" outflow boundary condition in practice may inherit some risk of ill-posedness. In Kracmar & Neustupa [59] this possible defect is treated by reformulating the problem as a variational inequality.
1.4. Discretization in space. In setting up a finite element model, one starts from the variational formulation of the problem: The choice of the function spaces H ⊂ H 1 (Ω) and L ⊂ L 2 (Ω) depends on the specific boundary conditions imposed in the problem to be solved. On a finite element mesh T h = {K} of Ω , one defines spaces of "discrete" trial and test functions, H h ⊂ H and L h ⊂ L. The "cells" K are (closed and convex) triangles or quadrilaterals in 2d and tetrahedra or hexahedra in 3d, and h ∈ R + is a mesh-size parameter, i.e., h := max K∈T h h K and h K := diam(K) . For a detailed description of the necessary structural regularity requirements on "admissible" meshes T h and mesh families {T h } h∈R+ , we refer to the standard literature, e.g., Ciarlet [30] and Brenner & Scott [26] (see Figure 2). Remark 4. In order to ease mesh refinement, we allow cells to have nodes which lie on midpoints of faces or edges of neighboring cells. But at most one such "hanging node" is permitted on each face or edge. There are no degrees of freedom associated to these irregular nodes as there the value of a finite element function is determined by interpolation of neighboring nodal values (see Carey & Oden [28] for more details).
The discrete analogues of (15), (16) then read as follows: where v in h is a suitable approximation of the inflow data v in , defined on all of Ω. The notation H h ⊂ H indicates that in this discretization the spaces H h may be "nonconforming", i.e., the discrete velocities v h are continuous across the interelement boundaries and zero along the rigid boundaries only in an approximate sense; in this case the discrete forms a h (·, ·) , b h (·, ·), n h (·, ·, ·) and the discrete "energy norm" ∇ · h are defined in the piecewise sense, In order that (17), (18) is a stable approximation to (15), (16), as h → 0, it is crucial that the spaces H h ×L h satisfy a compatibility condition, the so-called "infsup" or "Babuska-Brezzi" condition, which is the discrete analogue of the continuous inf-sup stability estimate (7): Here, the constant γ is required to be independent of h . This ensures that the problems (17) (0) The Q 1 /P 0 Stokes element For completeness, we mention the simplest (quadrilateral) Stokes element, the socalled Q 1 /P 0 element, which uses continuous bilinear velocity and piecewise constant pressure approximation; the triangular counterpart of this element is not consistent. This element pair is formally second-order accurate but lacks stability since the uniform inf-sup condition (19) is not satisfied on general meshes. Yet, this approximation has successfully been used in nonstationary flow computations (Gresho & Sani [44]). The explanation for the good performance of this "unstable" spatial discretization within a nonstationary computation is that the time-stepping, a pressure correction scheme, introduces sufficient stabilization as long as the time step is not too small.
(1) The nonconforming "rotated" d-linearQ 1 /P 0 Stokes element This is the natural quadrilateral analogue of the well-known triangular nonconforming finite element of Crouzeix/Raviart (Rannacher & Turek [79]). Its two-as well as three-dimensional versions have been implemented in state-of-the-art Navier-Stokes codes (Turek [97]). In two space dimensions, this nonconforming element uses piecewise "rotated" bi-linear reference shape functions for the velocities, spanned by {1, x, y, x 2 − y 2 } , and piecewise constant pressures.   (17), (18). We add certain least squares terms in the continuity equation (18) (so-called "pressure stabilization method"), where The correction terms on the right hand side have the effect that this modification is fully consistent, since the additional terms cancel out if the exact solution {v, p} of problem (15), (16) is inserted. In this way, one obtains a stable and consistent approximation of the Navier-Stokes problem (15), (16) with optimal-order accuracy (Hughes et al. [55], Becker & Rannacher [11]). However, from a practical point of view, this stabilization has several short-comings. First, being used together with the "free" outflow condition described above it induces a non-physical numerical boundary layer along the outflow boundary. Second, the evaluation of the various terms in the stabilization forms c h (χ h , p h ) and g h (v h ; χ h ) is very expensive, particularly in 3d. These problems can be resolved by using instead so-called "stabilization by pressure projection" or "local projection stabilization (LPS)" (Becker & Braack [6], Braack et al. [22]). Here, the stabilization forms are g h = 0 and where π 2h denotes the projection into a space L 2h defined on a coarser mesh T 2h . The resulting scheme is of second-order accurate, the evaluation of the system matrices is cheap and the consistency defect at the outflow boundary is avoided. Further, the resulting discretization has conceptional advantages in solving flow control problems by the indirect approach (based on the KKT system).
(3) Higher-order Stokes elements One of the main advantages of the Galerkin finite element method is that it provides systematic ways of dealing with complex geometries and of constructing higher-order approximation. Popular examples are: • The triangular P # 2 /P d 1 element with continuous quadratic velocity (augmented by a cubic bulb) and discontinuous linear pressures.
• The quadrilateral Q c 2 /Q c 2 element with continuous biquadratic velocity and pressure approximation and pressure stabilization.
• The triangular P c 2 /P c 1 and quadrilateral Q c 2 /Q c 1 Taylor-Hood elements with continuous quadratic/biquadratic velocity and linear/bilinear pressure approximation, or the even better element Q c 2 /P d 1 with discontinuous linear pressure approximation.
All these Stokes elements are inf-sup stable and third-order accurate for the velocity and second-order for the pressure, the errors measured in the respective L 2 norms. Practical experience shows that they yield much better approximations than the lower-order elements.

1.4.2.
Treating dominant transport. In the case of higher Reynolds number (e.g., Re > 500 for the 2d driven cavity, and Re > 50 for pipe-flow around a cylinder) the finite element models (17), (18) and (17), (20) may become unstable since they essentially use central-differences-like discretization of the advection term. This instability most frequently occurs in form of a drastic slow-down or even break-down of the iteration processes for solving the algebraic problems; in the extreme case the possibly existing "mathematical" solution contains nonphysical oscillatory components. In order to avoid these effects some additional stabilization is required. However, the use of simple first-order artificial viscosity or upwinding is not advisable since it introduces too much numerical diffusion.
The Galerkin structure of the FEM is fully respected by the so-called "streamline diffusion method". The idea of streamline diffusion is to introduce artificial diffusion acting only in the transport direction while maintaining the secondorder consistency of the scheme. This can be achieved either by augmenting the test space by direction-oriented terms resulting in a "Petrov-Galerkin method" with streamline-oriented upwinding ("SUPG method"), or by adding certain leastsquares terms to the discretization including streamline-oriented diffusion ("LSSD method", Hughes & Brooks [54], Franca & Frey [36], Tobiska & Verfürth [90], Braack et al. [22]).
In the following, we describe the full LSSD method for the Navier-Stokes system. To this end, we introduce the product Hilbert-spaces V := H × L of pairs u := {v, p} and ϕ = {ψ, χ} as well as their discrete analogues On these spaces, we define the semi-linear form and the linear functional F (ϕ) := (f, ψ). Then, the variational formulation (8) of the Navier-Stokes equations is written in the following compact form: For defining the stabilization, we introduce the matrix differential operator A(u) , depending on u = {v, p} , and the forcing vector F by Then, with the weighted L 2 scalar product the LSSD stabilized finite element approximation reads: This discretization contains several features. The stabilization form contains the terms where the first one stabilizes the pressure-velocity coupling for the Q 1 /Q 1 Stokes element, the second one stabilizes the transport operator, and the third one enhances mass conservation. The other terms introduced in the stabilization are correction terms which guarantee strong second-order consistency for the stabilized scheme. Practical experience and analysis suggest the following choice of the stabilization parameters: with values α ≈ 1 12 and β ≈ 1 6 . In practice, the terms −ν∆v h as well as −ν∆ψ h in the stabilization are usually dropped, since they vanish or almost vanish on the low-order elements considered. 1.4.3. The algebraic problems. The discrete Navier-Stokes problem including simultaneously pressure and streamline diffusion stabilization has to be converted into an algebraic system. To this end, we choose local "nodal bases" {ψ i h , i = 1, ..., N v } of the "velocity space" H h , and {χ i h , i = 1, ..., N p } of the "pressure space" L h , and expand the unknown x j ψ j h and p h = Np j=1 y j χ j h . Accordingly, we introduce the following matrices and vectors: . Notice that the nonhomogeneous inflow data v h|Γin = v in h is implicitly incorporated into the system, and that the stabilization only acts on velocity basis functions corresponding to interior nodes. Further, we introduce the velocity and pressure "mass matrices": With this notation the variational problem (22), can equivalently be written in form of an algebraic system for the vectors x ∈ R Nv and y ∈ R Np of expansion coefficients: This system has essentially the features of a saddle-point problem (since C small) and is generically nonsymmetric. This requires special care in designing efficient and robust iterative solvers such as PCG or multigrid methods (Vanka [98] and Turek [97]).

Discretization in time.
We now consider the nonstationary Navier-Stokes problem: Find v ∈ v in + H and p ∈ L , such that v(0) = v 0 and, for t > 0, where it is implicitly assumed that v as a function of time is continuous on [0, ∞) and continuously differentiable on (0, ∞). The choice of the function spaces H ⊂ H 1 (Ω) d and L ⊂ L 2 (Ω) depends again on the specific boundary conditions chosen for the problem to be solved. Due to the incompressibility constraint the nonstationary Navier-Stokes system has the character of a "differential-algebraic equation" (in short "DAE") of "index" two, in the language of ODE theory. This requires an implicit treatment of the pressure within the time-stepping process, while the other flow quantities may, in principle, be treated explicitly.
1.5.1. The Rothe Method. In the "Rothe Method", at first, the time variable is discretized by one of the common time differencing schemes. For example, the backward Euler scheme leads to a sequence of stationary Navier-Stokes-like problems for all {ϕ, χ} ∈ H × L , where k n = t n − t n−1 is the time step. Each of these problems is then solved by some spatial discretization method as described in the preceding section. This provides the flexibility to vary the spatial discretization, i.e. the mesh or the type of trial functions in the finite element method, during the time stepping process. In the classical Rothe method the time discretization scheme is kept fixed and only the size of the time step may change. The transfer between different spatial meshes is accomplished by some interpolation or projection method. Within the Rothe approach the most natural simultaneous discretization in space and time is the space-time Galerkin method, which is based on a corresponding space-time variational formulation of the problem, for all ϕ = (ψ, χ) T ∈ X , with the semi-linear form A(u)(ϕ) as defined above. The function space X is given as

ROLF RANNACHER
To introduce the semi-discretization in time, we partition the time intervalĪ The discretization parameter k is given as a piecewise constant function by k| Im = k m , for m = 1, . . . , M . On the subintervals I m , we define the following semi-discrete spaces X where P r (I m , Y ) denotes the space of polynomials up to degree r on I m with values in Y . To account for the possible discontinuity of a function u k at time points t m , we introduce the notation i. e., u + k,m and u − k,m are the limits "from above" and "from below" at time t m , respectively, while [u k ] m is the corresponding "jump" of u k at t m .
Then, the dG(r) ("discontinuous Galerkin") semi-discretization of problem (29) Remark 5. Due to the discontinuity of the test functions, the dG(r) discretization decouples into a time stepping scheme. For example, the dG(0) discretization is a variant of the backward Euler method, while the dG(1) discretization, after applying quadrature rules to the temporal integrals, corresponds to an implicit Runge-Kutta method.
The time-discrete problem (30) needs further discretization in space. For this, we can use any one of the spatial finite element spaces discussed above employing piecewise polynomial functions of degree s ≥ 1 . This spatial discretization is then referred to as a "(continuous Galerkin) cG(s) method" and the whole spacetime discretization as a "cG(s)dG(r) method" (Eriksson et al. [33]). To obtain the formulation of the fully discrete problem, we allow dynamic mesh change in time, but the local time steps k m are kept constant in space. To this end, with each time point t m , we associate a mesh T m h and corresponding (spatial) finite element spaces H sv,m h and L sp,m h . Then, we define the following space-time finite element space: Then, the fully discrete cG(s)dG(r) formulation of problem (29) This discretization can be used for inf-sup stable Stokes elements such as the Taylor-Hood element. Otherwise pressure stabilization has to be employed as described above. A similar modification is necessary in case of dominant transport. This results in the following modified fully discrete formulation: where The cellwise stabilization parameters α K,m and δ K,m are defined as with some constants α 0 and δ 0 . For details on the choice of these parameters, we refer to Franca & Frey [36] or Braack et al. [22]. In the computations below α 0 = δ 0 = 0.3.

1.5.2.
The Method of Lines. The more traditional approach to solving time-dependent problems is the "Method of Lines". At first, the spatial variable is discrete, e.g. by a finite element method as described above, leading to a DAE system of the form M 0 for t ≥ 0, with the initial value x(0) = x 0 . This DAE system is now discretized with respect to time. Let k denote the time-step size. Frequently used schemes belong to the simple "One-step-θ schemes" (step t n−1 → t n ): where x n ≈ x(t n ) and A n := A(x n ). Special cases are the "forward Euler scheme" for θ = 0 (first-order explicit), the most popular "Crank-Nicolson scheme" for θ = 1/2 (second-order implicit, A-stable), and the the "backward Euler scheme" for θ = 1 (first-order implicit, strongly A-stable). The most robust implicit Euler scheme is very dissipative and therefore not suitable for the time-accurate computation of nonstationary flow. In contrast to that, the Crank-Nicolson scheme has only very little dissipation but occasionally suffers from instabilities caused by the possible occurrence of rough perturbations in the data which are not damped out due to the only weak stability properties of this scheme (not strongly A-stable). This defect can in principle be cured by an adaptive step size selection but this may enforce the use of an unreasonably small time step.
Alternative schemes of higher order are based on the (diagonally) implicit Runge-Kutta formulas or the backward differencing multi-step formulas, both being well known from the ODE literature. These schemes, however, have not yet found wide applications in practical flow computations.
Another alternative is the so-called "Fractional-step-θ scheme" , and β = 1−α, in order to ensure second-order accuracy, and strong A-stability. Being strongly A-stable, for any choice of α ∈ (1/2, 1] , it possesses the full smoothing property in the case of rough initial data, in contrast to the Crank-Nicolson scheme (case α = 1/2). This has been shown for certain parameter ranges in Müller-Urbaniak [66]. Recently, a reduced variant of this scheme has been proposed and computationally analyzed in Glowinski et al. [43] (see also Glowinski [42]): Remark 6. The fractional-step-θ scheme was originally introduced as an operator splitting in order to separate the two main difficulties in solving problem (25) namely the nonlinearity causing nonsymmetry and the incompressibility constraint causing indefiniteness. In the first and last step one solves linear Stokes problems treating the nonlinearity explicitly, while in the middle step a nonlinear Burgers-type problem (without incompressibility constraint) is solved. The symmetric form of such a splitting scheme guarantees second-order accuracy of the approximation. However, nowadays, the efficient solution of the nonlinear incompressible Navier-Stokes equations is standard by the use of new multigrid techniques. Hence, the splitting of nonlinearity and incompressibility is no longer an important issue. Time-stepping schemes for the nonstationary Navier-Stokes equations have been computationally compared in Müller et al. [67] (see also Turek [95,97]).

Remark 7.
Another family of methods for time discretization are the so-called "projection methods". The classical example is the Chorin projection scheme (Chorin [29]), which has originally been designed in order to overcome the problem with the incompressibility constraint ∇ · v = 0 . Here, the continuity equation is decoupled from the momentum equation through an iterative process (again "operator splitting"). There are various schemes of this kind in the literature referred to as "projection method", "quasi-compressibility method", "SIMPLE method", etc.
For simplicity consider the case of pure homogeneous Dirichlet boundary conditions, v| ∂Ω = 0. Then, the simplest projection method reads as follows: For an admissible initial value v 0 , solve for n ≥ 1: (i) Implicit "Burgers step" forṽ n ∈ H: (ii) "Projection step" for v n :=ṽ n + k∇p n , wherep n ∈ H 1 (Ω) is determined by This time stepping scheme can be combined with any spatial discretization method, e.g., the finite element methods described above. Equation (35) amounts to a Poisson equation forp n with homogeneous Neumann boundary conditions. It is this non-physical boundary condition, ∂ np n | ∂Ω = 0, which has caused a lot of controversial discussion about the value of the projection method. Nevertheless, the method has proven to work well for representing the velocity field in many flow problems of physical interest. It is very economical as it requires in each time step only the solution of a (nonlinear) advection-diffusion system for v n (of Burgers equation type) and a scalar Neumann problem forp n . A rigorous convergence analysis shows that the quantitiesp n are indeed reasonable approximations to the pressure p(t n ) (Shen [86,87], Rannacher [72], E & Liu [32], Prohl [70]). The projection approach can be extended to formally higher order (Van Kan [99]) and is useful in constructing iterative solvers for the algebraic problems (Turek [96]).
1.6. Convergence analysis and error estimates. The usual error analysis for discretizations of the Navier-Stokes equations proceeds along the line of argument used for standard stationary or nonstationary nonlinear problems (see Ciarlet [30] for elliptic and Thomée [89] for parabolic problems). However, the nonstationary Navier-Stokes equations are a differential-algebraic system with all its inherent difficulties. The most essential of these difficulties origins from the loss of regularity at t = 0 due to the incompatibility of the initial data discussed in Section 1.2.1. The most complete a priori error analysis under natural regularity conditions on the solutions, u 0 ∈ H∩H 2 (Ω), has been given in the series of papers Heywood & Rannacher [48,49,51,52]. For second-order spatial (e.g., Q 1 /Q 1 elements) and first-order time discretization (e.g., backward Euler scheme) there holds the error estimate A deeper analysis exploiting the smoothing property of higher-order discretizations, yields smoothing error estimates of the form Since the loss of regularity at t = 0 is generic for the nonstationary Navier-Stokes equations, any error analysis of a higher-order discretization in space as well as in time has to be based on the smoothing property of the method. However, this critical aspect, which significantly complicates the analysis, is not often taken into account in the literature (Luskin & Rannacher [61,62]).
By counterexamples (Johnson et al. [57]) it is known that the smoothing property of standard numerical methods applied to nonlinear initial value problems is limited to some maximum order depending on the degree of nonlinearity and expected regularity of the solution at t = 0 . In the case of the Navier-Stokes equations this maximum-oder smoothing property is suspected to allow for the estimate for correspondingly higher-order spatial and time discretization. However, a detailed analysis of this conjecture is still missing. 2.1. Goal-oriented adaptivity. In the following, we discuss the principles of "goal-oriented" adaptivity in the Galerkin finite element method. Let the goal of a numerical simulation be the determination of a scalar quantity J(u) of a continuous model from its approximation by a discrete model, posed on a mesh T h = {K} consisting of small "cells". Then, the goal of adaptivity is to optimize the mesh T h (and model A h ) on the basis of an a posteriori error representation of the form Here, ρ K (u h ) are so-called "cell residuals" and ω K corresponding "cell weights", which describe the dependence of the error in the target quantity on variations of the cell residuals. These sensitivity factors are determined from the solution of certain "dual problems", as will be described in detail, below. Then, the mesh adaptation is oriented at the equilibration of the cell error indicators η K := ρ K (u h )ω K : In this context mesh adaptation for a set of physical quantities u = (u i ) n i=1 is based on "smoothness" or "residual" information like where D 2 h u i h is a second-order difference quotient and R i (u h ) a local residual, respectively. The complex error interaction has to be detected from a computed solution of a "dual" (or "adjoint") problem: i) the global error propagation in space and ii) the local interaction of physical error sources. On the bases of such weighted a posteriori error estimates, we will develop a feed-back process for generating economical meshes for the computation of the target quantity J(u) . This is the essence of the Dual Weighted Residual (DWR) Method.

2.2.
The theoretical framework. We develop the theoretical basis of the DWR method for goal-oriented a posteriori error estimation and mesh adaptation within an abstract setting. Let X be a function Hilbert space and L(·) a differentiable functional on X . We want to determine stationary points x ∈ X of L(·) by a Galerkin scheme in finite dimensional subspaces X h ⊂ X: Proposition 1 (Becker & Rannacher [14]). There holds the error representation with arbitrary y h ∈ X h , where the remainder R Proof. By elementary calculus using the trapezoidal rule with its integral reminder and Galerkin orthogonality, we conclude h , for arbirary y h ∈ X h . This proves the assertion.
For applying this general result to the Galerkin approximation in flow simulation, we start from the variational problem posed in a function space V , and its discrete Galerkin analogue in V h ⊂ V , Again the quantity of interest is given in terms of a functional J(·) . We introduce the Lagrangian functional L(·, ·) : V × V → R with "dual" variable z by satisfy the system The corresponding Galerkin approximations with the "primal" and "dual" residuals To embed this situation into the general framework underlying Proposition 1, we introduce the product spaces X := V × V, X h := V h × V h , with elements x := {u, z} ∈ X and x h := {u h , z h } ∈ X h , respectively, and the associated functional Then, at a stationary point x = {u, z} of L(·) there holds Using this relation, from Proposition 1, we can infer the following result.
Proposition 2 (Becker & Rannacher [14]). There hold the error representations with arbitrary ϕ h , ψ h ∈ V h , where the remainders R (2) h and R h are quadratic and cubic, respectively, in the errors e u := u−u h and e z := z−z h .

Application in transient simulations.
The efficient computation of instationary flows requires adaptivity in space as well as in time for target functionals of the kind J(u) = T 0 J 1 (u(t)) dt + J 2 (u(T )) . This adaptivity is based on a posteriori error estimates in which the residuals of spatial and time discretization are separated in order to allow for independent control of these partial discretizations, Here, space as well as time discretization is by Galerkin finite element schemes with continuous or discontinuous trial and test functions ("cG(r)" or "dG(r)" methods). For details of the discretization and the a posteriori error analysis, we refer to Eriksson et al. [33] and Schmich & Vexler [85] for standard parabolic problems, and to Besier & Rannacher [18] for the particular situation of the nonstationary Navier-Stokes equations. For deriving the a posteriori error estimates the error is split like where u k is the time-semi-discrete solution. The result of the basic Proposition 2 can be separately applied to both error components, i.e. at first to the time and then to the space discretization, which results in an error estimate of the form (51). This approach involves several steps of approximation: neglecting the higher-order remainder terms, replacing the arguments u k , z k in the time-residuals ρ(u k )(z − z k ) and ρ * (u k , z k )(u − u k ) by its spatially discrete analogues u kh , z kh , respectively, and approximating the continuous solutions u, z by local post-processing. For the algorithmic details and computational results, we refer to Besier & Rannacher [18]. For nonlinear problems the time-space adjoint problem involves the corresponding primal solution in its coefficients. This may cause serious storage problems particularly in long-time calculations on fine meshes, which requires the use of "data compression" or "check-pointing" for storage reduction (Becker et al. [7]).

2.4.
The main steps in the DWR method. In the following, we collect the main steps in the DWR method for a posteriori error estimation and mesh adaptation: -Solution of the discrete linear dual problem: usually corresponding to the work load of one additional Newton step. -Linearization by neglecting the high-order remainder terms R h , with a posteriori control if necessary. -Approximation of weights by local post-processing of the computed discrete solutions, e.g, patch-wise higher-order interpolation: [2]), with costs of a matrix-vector multiplication.
-Mesh adaptation on the basis of cell-wise error indicators η K by one of the common strategies: "error balancing", "fixed fraction", "fixed error reduction", or "look-ahead mesh optimization (see, e.g., Bangerth & Rannacher [2] or Rannacher [77]), with the goal of reaching most efficiently equilibration of the cell error indicators over the current meshes T h and eventually the desired error tolerance.

2.5.
Applications. The realization of the abstract error analysis developed above for the situation of the stationary Navier-Stokes problem leads to the following concrete form of the residuals (neglecting the contributions from pressure and transport stabilization): with the cellwise equation and normal-jump residuals (Becker [4]) This representation is directly evaluated and used in the stopping criterion of the adaptation process. For driving the mesh adaptation the following (non-negative) local indicators are used: In these identities the weights z v −z v h , z p −z p h and ∂ n z v are approximated by local higher-order post-processing as described in Subsection 2.4. The "dual" residual ρ * (z h )(u−u h ) and the corresponding local refinement indicators have an analogous form. The computational results presented below have been obtained on the basis of the error representation (49), which only involves the primal residual.
2.5.1. Drag computation in 2d stationary flow. For illustrating the features of the DWR method, we consider a 2-dimensional model situation involving the stationary Navier-Stokes equations for velocity/pressure (see Figure 4): The quantity of interest is the "drag coefficient" J(v, p) = c d defined by whereψ is an extension of ψ := (0, 1) T from S to Ω with support along S.
Here, σ(v, p) = ν(∇v + ∇v T ) + pI is the stress acting on the cylinder, D is its diameter, and U := max |v| is the reference inflow velocity. The Reynolds number is Re = U D/ν = 50, so that the flow is stationary. Notice that on the discrete level the two formulas (55) for c d differ. The second one is to be preferred as it yields significantly more accurate approximations. The discretization of this problem uses the equal-order Q 1 /Q 1 Stokes element on quadrilateral meshes. For mesh adaptation several different strategies are compared, which are driven by different local "error indicators" η K listed below: where the brackets [·] denote the jump of the enclosed quantity across intercell boundaries. The first three indicators (vorticity, pressure gradient, velocity gradient) merely measure "smoothness" of the discrete solution, the fourth one (residual) measures its consistency with the differential equation together with its smoothness (Oden et al. [68]), while the last one (weighted residual) additionally incorporates error sensitivity factors (Becker & Rannacher [12,13], Becker et al. [9], Machiels et al. [63]). The resulting adapted meshes shown in Figure 4 are quite different. Figure 5 shows the efficiency of the different refinement strategies. The weighted error indicator is clearly superior over the other approaches.

2.5.2.
Drag computation in 3d stationary flow. The second application concerns the computation of drag in a 3d configuration, namely channel flow around a cylinder with a square cross section, as shown in Figure 6 (Braack & Richter [24], John [56]). The model is chosen analogously to that in the prceding 2d example. For the discretization three different "Stokes elements" are considered: a) Q 2 /Q 1 -element -with uniform refinement, b) Q 1 /Q 1 -element -with local refinement driven by "smoothness" indicator, c) Q 1 /Q 1 -element -with local refinement driven by "weighted residual" indicator.
The corresponding results and adapted meshes are shown in Table 1 and in Figure 7, respectively. This example demonstrates the enormous potential of "model reduction" by sensitivity-driven mesh adaptation for achieving high accuracy at moderate cost.

2.5.3.
Drag computation in 2d nonstationary flow. As a third application, we consider a 2d analogue of the benchmark configuration "channel flow around a cylinder" but this time for a Reynolds number, which yields a nonstationary flow (Schäfer & Turek [84]). The configuration is shown in Figure 6. Now the target quantity is the temporal mean drag coefficient, For the further technical details of this test problem, we refer to Besier & Rannacher [18]. We employ the cG(1)dG(1) discretization in combination with the DWR method for simultaneous spatial and temporal mesh adaptation. The spatial meshes are kept constant over the whole time interval [0, T ] . The results of these computations are shown in Table 2, where we use the extrapolated value J d (u) = 1.6072872 as reference solution. We observe how the spatial and the temporal discretization errors are equilibrated and kept balanced under further refinement. Further, we notice a quite good agreement of the estimated and the actual discretization error especially on finer discretizations, i. e., the so-called "effectivity index" becomes Even though, in this example, we aim at efficiently computing the mean drag coefficient the maximum drag coefficient is also computed very accurately. Remark 8. We do not use dynamic spatial meshes in this example for two reasons. On the one hand, the use of dynamic meshes leads to wrong approximations of the drag coefficient if no additional projection steps are applied each time the spatial mesh is changed, which would be rather costly (Besier & Wollner [19]). The other reason becomes clear if we have a look at Figure 8 where the adaptive spatial mesh corresponding to the last line in Table 2 is shown. We observe that in order to precisely determine the mean drag coefficient it is not necessary to resolve the whole van Kármán vortex street. Only a small recirculation zone behind the obstacle is strongly refined. Since the vortices in this region develop relatively early, we may conclude that allowing dynamic meshes would not provide a notable reduction in the degrees of freedom needed to reach the same accuracy as with adaptively refined but fixed spatial meshes. In virtue of the additional effort on dynamic meshes due to more frequent matrix reassembling and the additional projection steps, we reach the conclusion that the use of dynamic spatial meshes does not make sense in this particular flow example.

2.6.
Balancing of discretization and iteration error. Usually the "exact" discrete solution u h is not available, but rather an approximationũ h obtained by some iterative method. This approximation spoils "Galerkin orthogonality", which is used in the derivation of the a posteriori error identity. Hence the goal is to design a strategy for controlling the (discrete) iteration based on a posteriori information, i.e., to balance the iteration error against the discretization error. In the context of duality-based error control such as underlying the DWR method, this problem has been analysed in a series of papers, in Becker et al. [16], Becker [3] for H 1 -and L 2 -error control in solving the Poisson and the Stokes problem, in Meidner et al. [65] for goal-oriented adaptivity, in Rannacher et al. [81] for eigenvalue problems, and in Rannacher & Vihharev [80] for nonlinear problems. For a survey of these results see Rannacher [77]. In the linear case, for discrete solutions u h satisfying Galerkin orthogonality, there holds For an only approximate discrete solutionũ h , which does not satisfies Galerkin orthogonality, we seek for an error representation of the form This question is answered by the following result (Meidner et al. [65]).
If a MG method is used on a sequence of sucessively refined meshes T 0 , T 1 , . . . , T L with canonical components, the following refined representation holds: Here,ẑ l ∈ V l can be chosen arbitrarily and R l (ṽ l ) are the iteration residuals on the mesh levels l = 1, . . . , L.
Proof. We only recall the argument for the first general error representation (58) and refer to Meidner et al. [65] for the proof of (59). There holds which already completes the proof.
Remark 9. We make the following remarks: -From the a posteriori error representation (58), we infer the following stopping criterion for the linear iteration (observe: ρ(u h )(ẑ h ) = 0 ): -The error representation (58) can be used for approximative solutionsũ h obtained by any solution process, where ρ(ũ L )(ẑ h ) may be considered as measuring the deviation ofũ h from satisfying Galerkin orthogonality. -The refined error representation (59) holds for V -, W -or F -cycles and for any type of smoothing. It allows not only to balance the iteration against the discretization error but also to tune the smoothing iteration separately on the different mesh levels.
2.6.1. Application in approximating the Stokes equations. The Stokes equation describes the behavior of a creeping incompressible fluid, where the notation from Subsection 1.1 is used. As model configuration, we consider the 2d channel flow around a circular obstacle as above where the quantity of interest is again the drag coefficient J(u) := c d . The discretization uses the equalorder (bilinear) Q 1 /Q 1 Stokes elements with pressure stabilization. The resulting discrete saddle point problems are solved by an MG method using the canonical mesh transfer operations and a "block ILU" smoothing (with 4+4 smoothing steps).
The computational results are shown in Tables 3, 4, and Figure 10. Here the total effectivity index is defined by We see that stopping the multigrid iteration in accordance to the discretization error reached on the current mesh (MG II) results in significant work savings compared to carrying the iteration down to round-off error level (MG I).   This has already been successfully demonstrated for various highly complex flow problems, e.g., including heat transfer and chemical reactions (Braack & Rannacher [23], Becker et al. [8]) and fluid-structure interaction (Dunne & Rannacher [31] and Bönisch et al. [20]). "Model adaptivity" has been considered in Braack & Ern [21].
3. Optimal flow control. In this third lecture, we discuss the use of the DWR method for adaptive discretization in optimal flow control and model calibration.
3.1. Thoughts on adaptivity in optimization. A general optimization problem seeks an optimal control q ∈ Q and an associated state u ∈ V , such that where J(·, ·) is a given cost functional defined on the state space V and control space Q . The notion of "admissibility" of states u = u(q) in PDE-constrained optimization requires some thoughts (Becker et al. [10]): -Discretization introduces perturbation in the state equation.
-Accuracy in discretization of PDEs is expensive.
-The efficient numerical solution of optimization problems governed by PDEs requires work reduction by adaptive discretization. -In PDEs the choice of error measures is a delicate matter: -Accuracy requirements should observe intrinsic problem sensitivities and is a modeling issue.
The goal is the development of a universal approach for dimension reduction by weakening admissibility requirements, which can be applied for general systems of PDEs and general (Galerkin) discretization and does not rely on mostly unknown (or worst-case oriented) coercivity properties.

3.2.
Theoretical framework (revisited). We want to apply the abstract result of Proposition 1 to the situation of functional minimization. With finite dimensional subspaces V h ×Q h ⊂ V × Q the Galerkin discretization of the optimization problem (63) seeks pairs {u h , q h } ∈ V h ×Q h , such that We introduce the corresponding Langrangian functional L(u, q, λ) := J(u, q) − A(u, q)(λ) with the adjoint variable λ ∈ V . Under certain conditions, which have to be guaranteed in the particular concrete situation considered, the solutions of the optimization problem (63) are given by stationary points of L . These satisfy the continuous optimality systems, the so-called Karush-Kuhn-Tucker (KKT) systems: The corresponding Galerkin approximation, defining stationary points of L on The associated residuals are given by For controlling the quality of this approximation the natural choice is to measure its error with respect to the cost functional J(·, ·) . In this case, the general Proposition 1 yields the following result.

Remark 10.
In approximating parameter estimation problems the choice of the least-squares cost functional J(·, ·) is questionable for a posteriori error control and mesh design. Consider, for example, the model problem where the goal is to determine the coefficient q ∈ R by comparing the resulting state u(q) with given measurementsū , J(u, q) := 1 2 u−ū 2 + 1 2 α q 2 → min, 0 < α 1.
The corresponding KKT system reads −(ϕ, u−ū) + (∇ϕ, ∇λ) + (qϕ, λ) = 0 ∀ϕ ∈ V, where V := H 1 0 (Ω) and Q := R . In this case the error identity based on the "artificial" least-squares cost functional J(u, q) may be useless for steering the mesh adaptation. Indeed, if the parameter q ≥ 0 is identifiable, we have u =ū at the stationary point, and the adjoint variable λ satisfies −∆λ + qλ = u−ū = 0, which implies that λ ≡ 0 . This renders the error representation useless for steering the mesh adaptation process. The difficulty can be handled in different ways: -An energy-norm-type a posteriori error estimate for |q − q h | can be derived based on coercivity estimates for the saddle-point problem. However, the stability constant in this estimate is usually unknown and depends on α −1 . -An a posteriori error estimate for a suitable expression of q−q h can be derived by an "outer" duality argument (Becker & Vexler [15]).

3.3.
Algorithmical issues. Next, we discuss some computational issues in solving optimization problems. a) The direct solution approach: We denote by S(·) : Q → V the solution operator of the state equation, A(S(q), q)(ψ) = 0 ∀ψ ∈ V, where the local existence and sufficient regularity of S(q) is assumed. Then, we reformulate (63) as unconstrained optimal control problem, j(q) := J(S(q), q) → min, with the "reduced cost functional" j(·) : Q → R . The first and second-order necessary conditions for an optimal q are It is usually computed by the so-called Newton-SQP method, with update q n+1 := q n + δq n . The single Newton steps are solved by an outer Krylov-space method with inner multigrid acceleration. This requires the evaluation of the directional derivatives j (q n )(δq n , χ) and j (q n )(χ) for fixed χ, which are given in terms of the Lagrangian functional by j (q)(δq) = L q (u, q, λ)(δq), j (q)(δq, δr) = L qq (u, q, λ)(δq, δr) + L uq (u, q, λ)(δu, δr) where for given δq, δr the quantities δu, δλ are obtained by solving certain auxiliary linear problems. The above relations are formulated on the continuous level. In practice they have to be carried out on the discrete level in appropriate finite element The indirect solution approach: The optimal control q ∈ Q and associated state u(q) ∈ V are simultaneously approximated by solving the fully coupled discretized KKT system (68).
This may be done similarly as above by an outer Krylov-space method with multigrid accelleration. There are various pros and cons of the indirect approach, particularly in nonstationary situations (Meidner [64]).
Remark 11. We remark on additional computational aspects and extensions.
-The iterative solution of the discretized KKT system (77)-(79) requires proper balancing of the discretization and iteration errors. This can again be achieved on the basis of corresponding goal-oriented a posteriori error estimates (Meidner et al. [65]). -The solution of nonstationary optimization problems together with goal-oriented mesh adaptation requires special care. In this case the primal state occurs in the coefficients of the adjoint problem, which leads to high storage requirements. For the treatment of this problem, we refer to Becker et al. [8] and Meidner [64].
-The incorporation of additional control and state constraints into the DWR method, such as lower bounds on the lift coefficient, c lift ≥ c 0 , or the suppression of local recirculation, v 1 ≥ 0 is possible. For first steps into this direction, we refer to Wollner [102] and Rannacher et al. [78]. -As a result of mesh adaptation by the DWR method the optimization cycle is performed on strongly nonuniform and partially rather coarse meshes. This may lead at the end to "optimal" controls q opt h and states u opt h , which are far from being "admissible" in the strict sense. If more admissibility is required, then an improved stateũ opt h may be recovered by approximating the state equation in a larger spaceṼ h ⊃ V h defined on a uniform finer mesh: Applications to flow problems.

3.4.1.
Flow control: drag minimization (Becker [5] and Becker et al. [7]). We again consider the configuration of the 2d flow benchmark "channel flow around a cylinder" described above. The goal is to minimize the drag coefficient by implying a pressure drop at two outlets Γ Q at the top and the bottom of the channel as shown in Figure 11. The Reynolds number of the uncontrolled flow is Re =Ū 2 D/ν = 40 such that initially stationary flow occurs. Therefore, the stationary Navier-Stokes equations are taken as mathematical model for the optimization process. The state equation reads in variational form as follows: where the governing semilinear form is defined by The results of the optimization are shown in Figure 12 and Table 5. We see that using the DWR approach for mesh adaptation the optimization takes place on meshes, which resolve the flow only in the vicinity of the obstacle (where the drag "lives"), the control boundaries, and the inflow boundary. It seems clear that this leads to significant work savings compared to a computation on a corresponding uniformly refined mesh. This is well confirmed by the data in Table 5. However, looking at the structure of the flow field corresponding to the optimal control in Figure 12, we must be concerned about the stability, i.e., physical relevance, of this stationary "optimal" state. This question will be further addressed by conducting a numerical stability analysis, below.   [75]). As second example, we consider a problem of "model calibration" for the stationary Navier-Stokes equations in case when the flow is driven by an unknown pressure drop imposed at parts of the boundary. This is a typical situation in modeling blood flow through human blood vessels. In Subsection 1.3, we have discussed that the choice of proper boundary conditions at "artificial" outflow boundaries is a delicate problem and actually part of the modelling process. In case that certain measurement data of the flow to be simulated are available such boundary conditions, i.e., mean pressure values, may be determined by numerical model calibration.
As an example, we consider the 2d configuration shown in Figure 13. Here, the correct Dirichlet or Neumann data at the two inlets Γ (i) in , i = 1, 2 , and the outlet Γ out are usually unknown. But this knowledge is necessary for determining the correct flow in the considered subdomain, which may be only part of a much larger system of blood vessels. For instance, the model is well-posed if the mean pressures q (i) at the inlets Γ (i) in , i = 1, 2 , are known, while fixing that at the outlet Γ out to be zero.
The pressure values q i are to be determined from given measurement data, which are, for example, mean fluxes along the lines S i ∈ Ω as indicated in Figure 13, The optimization problem then consists in minimizing the "misfit functional" on the set of solutions to the flow model (81). The results obtained by using the DWR approach in this parameter identification process are shown in Figures 13  and 14. The superiority of sensitivity-driven over merely smoothness-oriented mesh adaptation is clearly seen.  4. Numerical stability analysis. In this fourth lecture, we consider the numerical stability analysis of stationary flows employing the concepts of linearized stability (Heywood & Rannacher [50]) and pseudospectra (Trefethen & Embree [92]). In particular, we will discuss the efficient computation of critical eigenvalues of the (nonsymmetric) stability eigenvalue problem associated to the linearized Navier-Stokes operator and corresponding pseudospectra.
with the tangent form Then, the stability of this base solution is determined by the growth property, of the solution operator S(t) : J 0 → J 0 of (82) where J 0 is the space of L 2 velocity fields, which are solenoidal and satisfy the boundary condition n · v| Γ rigid ∪Γin = 0 in the distributional sense. We introduce the space J 1 = {v ∈ H, ∇ · v = 0} and then have J 0 = cl · J 1 . In the relation (83) λ is the most critical eigenvalue of the non-symmetric eigenvalue problem associated to the tangent form a (v)(·, ·) : Find u := {v, p} ∈ V and λ ∈ C , such that Clearly, if one of these eigenvalues has real part Reλ < 0 , there are arbitrarily small perturbations of the base solution, which initially grow exponentially rendering the base solution "unstable". On the other hand, even if all eigenvalues satisfy Reλ > 0 this does not necessarily imply that small perturbations remain bounded. This is due to the possibly large size of the amplification constant A in (83) caused by the non-normality of the linearized Navier-Stokes operator A := A (v) . The spectrum, i.e. the set of eigenvalues, of A is denoted by Σ(A) . The variational form of the eigenvalue problem (84) and its adjoint analogue read where m(ψ, ϕ) = (ψ v , ϕ v ) . If m(u, u * ) = 0 , the boundary value problem has a solutionũ ∈ V ("generalized eigenfunction") and the eigenvalue λ has "defect" α(λ) ≥ 1 .

4.1.2.
The effect of non-normality and pseudospectra. The generic scenario in numerical stability analysis is that there are r ≥ 1 discrete eigenvalues {λ (i) h , i = 1, ..., r} which (due to asymmetries caused by discretization) are usually all simple and (for h → 0 ) approximate some eigenvalue λ ∈ Σ(A), such that (i) λ has geometric multiplicity m and hence trivial defect, or (ii) λ has geometric multiplicity less than r and hence defect α(λ) ≥ 1 . The existence of an eigenvalue with 0 < Reλ 1 but α(λ) ≥ 1 may cause dynamic instability of the base flowû even for very small perturbations. This is reflected by the growth property with a certain rate, can be used as an indicator for α(λ) ≥ 1 . In this case the blowup behavior of the solution operator S(t) implies Consequently, for 0 < Reλ 1 the amplification constant A in the stability estimate becomes very large.
A similar destabilizing effect can also occur if λ is non-deficient. This is related to the concept of "pseudospectrum" of Trefethen et al. [93,94]. For ε ∈ R + the ε-pseudo-spectrum Σ ε ⊂ C of A (v) is defined by Remark 12. Pseudo-spectra are interesting only for non-normal matrices (closed operators in Banach spaces) since the pseudo-spectrum of a normal operator is the union of ε-circles around its eigenvalues.
The following result is related to the "easy part" of the so-called "Kreiss matrix theorem" (Trefethen & Embree [92] and Gerecht et al. [39]).
Hence a critical pseudo-spectrum may trigger instability even if all eigenvalues have positive real parts.

4.2.1.
Adaptive finite element discretization. The discretization of the stability eigenvalue problems (85) and (86) may again employ a stabilized finite element method (Q 1 /Q 1 Stokes element with GLS stabilization) using the bilinear form Then, the discrete primal and dual eigenvalue problems seek with the normalization m(u h , u h ) = m(u h , u * h ) = 1 . Remark 13. In computing potentially critical eigenvalues of A (v) , i.e. those with 0 < Reλ 1 , the following points have to be observed: -Reliable and efficient computation of critical eigenvalues requires i) error control by a posteriori error estimates, and ii) work reduction by "goal-oriented" mesh adaptation.
where a(·, ·) is an elliptic not necessarily symmetric sesquilinear form, which is compact w.r.t. the L 2 scalar product (·, ·). The corresponding Galerkin approximation in subspaces H h ⊂ H reads Then, under the consistency assumption abstract operator theory arguments yield the following a priori error estimates (Bramble & Osborn [25]): Corresponding a posteriori error estimates for eigenvalue problems of the type (94) have been derived in Heuveline & Rannacher [45], with the primal and dual cell residuals defined by This estimate is proved via "energy arguments". It usually crudely overestimates the true eigenvalue error since it does not correctly represent the interaction of primal and dual error sources, which is the essence of the DWR method described above. Using this more sophisticated approach, we obtain the following error representation: is cubic in the errors e h := u−u h and e * h := u * −u * h . The primal and dual residuals are defined by

4.2.3.
Approximation of the stability eigenvalue problem. Next, we consider the discrete primal and dual stability eigenvalue problems (92) and (93), where the approximate (stationary) base solutionû h is obtained from the discretized Navier-Stokes problem (17). For this setting, we have the following result.
Proposition 6 (Heuveline & Rannacher [46]). There holds the error representation for arbitraryψ h ,φ h , ψ h , ϕ h ∈ V h . The remainder R The proof uses the general result of Proposition 2. To this end, we embed the present situation into the framework of variational equations by introducing the product spaces Using this notation the base equation and the associated stability eigenvalue problem together with their discrete analogues can be written in compact form as The error in this galerkin approximation is measured in terms of the functional The corresponding continuous and discrete dual problems seek Z = {ẑ, z, π} ∈ V and A careful examination of these problems shows that {z, π} = {u * , λ} solves the corresponding adjoint eigenvalue problem, andẑ =û * is determined by This proves the asserted error representation.
The residual terms in the error representation (100) can again be split into the single cell distributions, analogously as in the representations (98) or (53), resulting in estimates of the form From this, we infer the following eigenvalue error estimate and a balancing criterium, which has to be satisfied in order to guarantee the reliability of the eigenvalue approximation: 4.2.4. Application to the drag minimization problem. We apply the strategies described above for analyzing the stability of the optimal stationary state obtained in Subsection 3.4.1. The resulting adapted meshes are shown in Figure 15. We see that the reliable eigenvalue computation requires more global mesh refinement than the optimization process. This reflects the fact that eigenvalues are closely related to the energy norm of eigenfunctions. The interplay between accuracy in computing the base solution and that in approximating the associated eigenvalues is depicted in Figure 16. We see that on coarse meshes the error in the base solution dominates but under mesh refinement based on (108) it falls below the eigenvalue error.  and that in approximating the associated eigenvalues ( * ).
Real and imaginary part of the computed critical eigenvalue depending on the prescribed control pressure are shown in Table 17. The optimal pressure is about q opt = 0.5 , which leads to negative Reλ < 0 , i.e., the corresponding optimal state v opt , turns out to be unstable as suspected.
since the "smallest" eigenvalues are to be computed. In the following, we summarize the ingredients of the algebraic solution process (Trefethen [91], Gerecht et al. [39]): -Arnoldi method using Krylov spaces K k = span{u 1 , Tu 1 , . . . , T k−1 u 1 } starting from a suitable vector u 1 (Saad [83]). -Computation of v j := T j u 1 , j = 1, . . . , k − 1 , by geometric MG on sequences of locally adapted meshes (most cost critical part). -Orthonormalization of Krylov vectors by stabilized Gram-Schmidt or Householder algorithm. -Reduction to k-dimensional Hessenberg matrix H k ∈ R k×k (k N ): -Approximate eigenvalues λ h,k are obtained by the k-dimensional eigenvalue problem H kṽh,k = λ −1 h,kṽ h,k , which can be solved by the QR-method in O(k 2 ) operations.
-Corresponding eigenvectors are obtained by a modified QR decomposition of the singular system and by setting v h,k := V kṽh,k .
A more detailed description of this algorithm is contained in Gerecht et al. [39]. For an alternative approach also employing multigrid techniques, we refer to Bertsch & Heuveline [17].
For stopping the Arnoldi iteration at the level of the discretization error, we use a criterion based on the following a posteriori result.
For reasonable iteration number k the differences may be rather large and therefore, using η it k := ρ(u h,k , λ h,k )(u * h,k ) as an indicator for the iteration error, requires a Then, the computation of the singular value decomposition of the Hessenberg matrix H h,k or its inverse H −1 h,k yields approximations to Σ ε (T h ) or Σ ε (A h ) , respectively. In turn, these approximate Σ ε (T ) or Σ ε (A (û)) , as h → 0 , where T = A (û) −1 . This is guaranteed by the following result.
Proposition 9 (Gerecht et al, [39]). (i) Suppose that all eigenvalues of the linearized Navier-Stokes operator have positive real parts, i.e., Σ(A (v)) ⊂ C + := {z ∈ C, Rez > 0} . Then, for sufficiently small h, also Σ(A h (v)) ⊂ C + . (ii) The pseudospectra of the discrete operators T h converge to that of the continuous operator T = A (û) −1 in the sense that, for sufficiently small h > 0 , (iii) The pseudospectra of the discrete operators A h converge to those of the continuous operator A = A (û) in the sense that, for ε(1 + c * h 2 ) ≤ 1, 4.4. Computational results.

4.4.1.
Pseudospectra of Burgers and Navier-Stokes operators. At first, we consider pseudospectra of the (vector) Burgers operator compared to those of the Navier-Stokes operator both linearized about Couette flow as base flow. The results for the case of Neumann-type inflow boundary conditions are shown in Figure 18. We see a significant difference in the distribution of eigenvalues and structure of pseudospectra. Next, we want to test the relevance of the "deficiency test" developed from the result of Proposition 8. For this, we choose ε = 2Reλ crit v * −1 1 and check whether with the computed pseudospectrum of A (v), there holds λ ε := λ crit − ε v * = −λ crit ∈ Σ ε (A (v)).
(120) Table 6 contains the critical eigenvalue λ crit and the "test quantity" v * h depending on the Reynolds number. For Re = 350 the critical eigenvalue is λ (350) crit ≈ 0.0282 . For v * ≈ 450 and ε = 2Reλ crit v * −1 ≈ 1.3 · 10 −4 , there holds 4.4.2. Pseudospectra in the channel flow benchmark. Next, we consider the pseudospectra of the Navier-Stokes operator in the configuration of the flow benchmark problem "channel flow around a cylinder" linearized about its stationary solution. The configuration and stationary base flow are depicted in Figure 20. The viscosity is chosen such that the Reynolds number is small enough, 20 ≤ Re ≤ 40 , to guarantee stationarity of the base flow. Already for Re = 60 the flow is nonstationary (time periodic). The stationary base flow is computed using Dirichlet inflow conditions but the associated eigenvalue problem of the linearized Navier-Stokes operator is considered with Neumann ("free") inflow conditions. Table 7 contains the results. In the case of perturbations required to satisfy Dirichlet inflow conditions the stationary base flow is stable up to Re = 45 . In the present case of perturbations without any inflow constraints (free inflow) at Re = 40 the critical eigenvalue has positive but very small real part, Reλ min ≈ 0.003. Hence, the precise stability analysis requires the determination of the corresponding pseudospectrum. This is shown in Figure 21. Though, for Re = 40 the real part of the most critical (positive) eigenvalue is rather small, the corresponding 10 −2 -pseudospectrum reaches only a little into the negative complex half-plane. This is in good agreement with the development of the quantity v * h under mesh refinement shown in Table 8. Table 7. Eigenvalues with smallest real part of the linearized ("channel flow") Navier-Stokes operator with "free inflow" conditions for different Reynolds numbers and increasing level of refinement.