Error estimation and adaptive mesh refinement for aerodynamic flows Ralf Hartmann1 and Paul Houston2 1 Institute of Aerodynamics and Flow Technology DLR (German Aerospace Center) Lilienthalplatz 7, 38108 Braunschweig,Germany [email protected] 2 School of Mathematical Sciences University of Nottingham University Park, Nottingham, NG7 2RD, UK [email protected] Contents 1 Introduction 4 1.1 Elements of function space theory. . . . . . . . . . . . . . . . . . . . . . . . . 6 1.1.1 Spaces of continuous functions . . . . . . . . . . . . . . . . . . . . . . 6 1.1.2 Spaces of integrable functions . . . . . . . . . . . . . . . . . . . . . . . 7 1.1.3 Sobolev spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2 Motivation: Linear problems and adjoint equations 9 2.1 Error estimation for linear problems . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Derivation of adjoint problems for linear primal problems . . . . . . . . . . . 11 2.3 The linear advection equation and adjoint problems . . . . . . . . . . . . . . 11 2.4 Numerical example: Linear advection equation . . . . . . . . . . . . . . . . . 13 3 Discontinuous Galerkin methods for compressible flows and their corre- sponding adjoint problems 15 3.1 The compressible Euler equations . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 Derivation of adjoint problems for nonlinear primal problems . . . . . . . . . 17 3.3 The adjoint equations to the compressible Euler equations . . . . . . . . . . . 17 3.4 DG discretization of the compressible Euler equations . . . . . . . . . . . . . 18 3.5 Consistency and adjoint consistency . . . . . . . . . . . . . . . . . . . . . . . 21 3.6 The compressible Navier-Stokes equations . . . . . . . . . . . . . . . . . . . . 22 3.7 The adjoint equations to the compressible Navier-Stokes equations . . . . . . 23 3.8 DG discretization of the compressible Navier-Stokes equations . . . . . . . . . 24 3.9 Consistency and adjoint consistency . . . . . . . . . . . . . . . . . . . . . . . 26 1 4 Adjoint-based error estimation and adaptive mesh refinement 28 4.1 Error estimation and mesh refinement for single target quantities . . . . . . . 28 4.2 Error estimation for multiple target quantities . . . . . . . . . . . . . . . . . . 31 4.2.1 The standard approach . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2.2 A new approach . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 4.3 Adaptive refinement for multiple target quantities . . . . . . . . . . . . . . . 33 4.4 Derivation of residual-based indicators . . . . . . . . . . . . . . . . . . . . . . 35 4.5 Numerical examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 4.5.1 Ringleb flow problem. . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.5.2 Supersonic flow past a wedge . . . . . . . . . . . . . . . . . . . . . . . 37 4.5.3 Supersonic flow past a BAC3-11 airfoil . . . . . . . . . . . . . . . . . . 41 4.5.4 Supersonic viscous flow around the NACA0012 airfoil . . . . . . . . . 45 4.5.5 Comparison of the approximate error representation for viscous and inviscid flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 4.5.6 Errorestimationandadjoint-basedrefinementformultipletargetquan- tities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 5 Development of anisotropic mesh adaptation 63 5.1 Model problem and discretization . . . . . . . . . . . . . . . . . . . . . . . . . 64 5.2 Meshes, finite element spaces and traces . . . . . . . . . . . . . . . . . . . . . 65 5.3 Interior penalty discontinuous Galerkin method . . . . . . . . . . . . . . . . . 66 5.4 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.5 Approximation results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 5.6 A priori error bounds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 5.7 A posteriori error estimation and adaptivity . . . . . . . . . . . . . . . . . . . 81 5.8 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.8.1 Singularly perturbed advection-diffusion problem . . . . . . . . . . . . 83 5.8.2 ADIGMA MTC3: Laminar flow around a NACA0012 airfoil . . . . . . 87 5.8.3 ADIGMA BTC0: Laminar flow around streamlined body . . . . . . . 90 6 High-order/hp–adaptive finite element methods for compressible flows 91 6.1 Model problem and discretization . . . . . . . . . . . . . . . . . . . . . . . . . 93 6.1.1 Meshes and finite element spaces . . . . . . . . . . . . . . . . . . . . . 93 6.1.2 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 6.2 hp-Error bounds on the hypercube . . . . . . . . . . . . . . . . . . . . . . . . 96 6.2.1 Isotropic polynomials degrees . . . . . . . . . . . . . . . . . . . . . . . 97 6.2.2 Anisotropic polynomial degrees . . . . . . . . . . . . . . . . . . . . . . 98 6.3 A priori error analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 6.4 hp–Adaptivity on isotropically refined meshes . . . . . . . . . . . . . . . . . . 106 6.4.1 hp–extension control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 6.5 Numerical Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 6.5.1 Mixed hyperbolic–elliptic problem . . . . . . . . . . . . . . . . . . . . 111 6.5.2 ADIGMA MTC1: Inviscid flow around a NACA0012 airfoil . . . . . . 117 6.5.3 ADIGMA MTC3: Laminar flow around a NACA0012 airfoil . . . . . . 120 6.6 Anisotropic hp–mesh adaptation . . . . . . . . . . . . . . . . . . . . . . . . . 121 6.7 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 2 6.7.1 Singularly perturbed advection–diffusion problem . . . . . . . . . . . . 124 6.7.2 Mixed hyperbolic–elliptic problem . . . . . . . . . . . . . . . . . . . . 126 6.7.3 ADIGMA MTC3: Laminar flow around a NACA0012 airfoil . . . . . . 129 7 Application of error estimation and adaptation to complex flows 132 7.1 ADIGMA BTC0: Laminar flow around streamlined body . . . . . . . . . . . 133 7.2 ADIGMA BTC3: Laminar flow around delta wing . . . . . . . . . . . . . . . 137 7.3 ADIGMA BTC1: L1T2 high-lift configuration . . . . . . . . . . . . . . . . . . 142 7.4 ADIGMA BTC0: Turbulent flow around streamlined body . . . . . . . . . . 149 7.5 ADIGMA CTC4 (modified): Subsonic turbulent flow around DLR-F6 wing- body configuration without fairing . . . . . . . . . . . . . . . . . . . . . . . . 150 Acknowledgements 153 Bibliography 157 3 1 Introduction Computational fluid dynamics (CFD) has become a key technology in the development of new products in the aeronautical industry. During the last decade aerodynamic design engi- neers have progressively adapted their way-of-working to take advantage of the possibilities offered by new CFD capabilities based on the solution of the Euler and Navier–Stokes equa- tions. Significant improvements in physical modelling and solution algorithms have been as important as the enormous increase of computer power to enable numerical simulations at all stages of aircraft development. However, despite the progress made in CFD, in terms of user time and computational resources, large aerodynamic simulations of viscous flows around complex configurations are still very expensive. The requirement to reliably compute results with a sufficient level of accuracy within shortturn-aroundtimes places severe constraints on the application of CFD. Indeed, within CFD the most popularclass of methods which are currently used in industrial codes are based on employing finite volume methods. While in principal these methods are second–order accurate, in practice their convergence order deteriorates to somewhere between first– and second–order on irregular and/or highly stretched meshes. Thereby, for reliablenumericalpredictionstobemadebysuchmethods,extremelyfinemesheswithalarge numberof degrees of freedomarerequired,which inturnleads to excessively largecomputing times. As an alternative approach, in recent years there has been significant interest in the development of high–order discretization methods; this is particularly evidenced by the funding of the EU Framework 6 project ADIGMA [82] (Adaptive higher order variational methods for aerospace applications) comprising of a consortium of academic and industrial partners. On a given mesh they allow for an improved prediction of critical flow phenomena, such as boundary layers, wakes, and vortices, for example, as well as force coefficients, e.g., drag, lift, moment. In particular, high–order methods are capable of achieving the same level of accuracy while exploiting significantly fewer degrees of freedom compared with classical finite volume methods. One extremely promising class of high-order schemes based on the finite element frame- work are Discontinuous Galerkin (DG, for short) methods. Indeed, the development of DG methods for the numerical approximation of the Euler and Navier-Stokes equations is an extremely exciting research topic which is currently being developed by a number of groups all over the world, cf. [14, 15, 19, 20, 34, 38, 39, 50, 59, 61, 62, 95, 107, 108], for example. DG methodshave several importantadvantages over well established finitevolume methods. The concept of higher-order discretization is inherent to the DG method. The stencil is minimal in the sense that each element communicates only with its direct neighbors. In particular, in contrast to the increasing stencil size needed to increase the accuracy of classical finite volume methods, the stencil of DG methods is the same for any order of accuracy which has important advantages for the implementation of boundaryconditions and for the parallel efficiency of the method. Moreover, due this simple communication at element interfaces, elements with so–called hanging nodes can be easily treated, a fact that simplifies local mesh refinement(h–refinement). Additionally, thecommunication atelementinterfaces is identical for any order of the method which simplifies the use of methods with different polynomial orders p in adjacent elements. This allows for the variation of the order of polynomials over the computational domain (p–refinement), which in combination with h–refinement leads to so–called hp–adaptivity. 4 Mesh adaptation in finite element discretizations should be based on rigorous a posteriori error estimates; for hyperbolic/nearly–hyperbolic equations such estimates should reflect the inherentmechanismsoferrorpropagation(see[70,76]). Theseconsiderationsareparticularly important when local quantities such as point values, local averages or flux integrals of the analytical solution are to be computed with high accuracy. In the context of aerodynamic flow simulations, it is of vital importance that certain force coefficients, such as the drag, lift and moment on a body immersed within a compressible fluid, are reliably and efficiently computed. Selective error estimates of this kind can be obtained by the optimal control technique proposed in [36] and [23] which is based on duality arguments analogous to those from the a priori error analysis of finite element methods. In the resulting a posteriori error estimates the element-residuals of the computed solution are multiplied by local weights involving the adjoint solution. These weights represent the sensitivity of the relevant error quantity withrespecttovariations of thelocal meshsize. Sincetheadjoint solution is usually unknown analytically, it has to be approximated numerically. On the basis of the resulting a posteriori error estimate the current mesh is locally adapted and then new approximations to the primal and adjoint solution are computed. This feed-back process is repeated, for instance, until the required error tolerance is reached. In this way, optimal meshes, or in the hp–setting, optimal finite element spaces can be obtained for various kinds of error measures, where optimal can mean most economical for achieving a prescribed accuracy TOL or most accurate for a given maximum number N of degrees of freedom. This approach is quite max universal as it can, in principle, be applied to almost any problem, linear or nonlinear, as long as it is posed in a variational setting. This lecture course covers the theory of so-called duality-based a posteriori error esti- mation of DG finite element methods. In particular, we formulate consistent and adjoint consistent DG methods for the numerical approximation of both the compressible Euler and Navier–Stokes equations; inthelattercase,theviscoustermsarediscretizedbasedonemploy- ing an interior penalty method. By exploiting a duality argument, adjoint–based a posteriori errorindicators willbeestablished. Moreover, application ofthesecomputableboundswithin automatic adaptive finite element algorithms will be developed. Here, a variety of isotropic and anisotropic adaptive strategies, as well as hp–mesh refinement will be investigated. The outline of these notes is as follows. In Section 2 we give an introduction to the adjoint–based a posteriori error estimation and mesh refinement for linear problems, and their subsequent exploitation within an automatic adaptive finite element algorithms. Then, in Section 3 we introduce both the compressible Euler and Navier–Stokes equations and formulate DG numerical methods for their discretization. In particular, here we will be con- cerned with the derivation of so-called adjoint consistent methods, which ensure the optimal approximation of target functionals of the underlying solution. Section 4 is devoted to the derivation of adjoint–based a posteriori error boundsfor the computed error in a given target functional of interest. Moreover, extensions to the case when there are multiple quantities of interest will be considered. The practical performance of these a posteriori error estimates within adaptive finite element algorithms will be studied through a series of numerical ex- periments. In Section 5 we consider the generalization of the above ideas to the case when anisotropicmeshrefinementis permitted. Inthissetting, wederivebotha priori anda poste- riori errorboundsfor theDGapproximation oflinear functionalsoftheunderlyinganalytical solution. The a priori analysis is fully explicit in terms of the anisotropy of the underlying computational mesh. Further, we introduce an anisotropic refinement algorithm, based on 5 choosing the most competitive subdivision of a given element from a series of trial (Carte- sian) refinements. The extension of these ideas to general anisotropic hp–version DG finite element methods is undertaken in Section 6. Finally, Section 7 is devoted to the application of goal-oriented adaptive finite element algorithms to complex aerodynamic flows, including three dimensional laminar flows as well as two and three dimensional turbulent flows. Before we embark on Section 2, we first take a brief excursion into the theory of function spaces to introduce the notational conventions used throughout these lecture notes. 1.1 Elements of function space theory The aim of this section is to provide a brief overview of some elementary results from the theory of function spaces and to introduce the notation which will be used throughout. For proofs and further technical details on classical function spaces the reader is referred to the monograph of Adams [1], for example. 1.1.1 Spaces of continuous functions Let N denote the set of all nonnegative integers. An n–tuple α = (α ,...,α ) in Nn will be 1 n referred to as a multi-index; the nonnegative integer α = α + + α is called the length 1 n | | | | ··· | | of the multi-index α. We define ∂α = ∂α1...∂αn, where ∂ = ∂/∂x for j = 1,...,n. 1 n j j Suppose that ω is an open set in Rn. For k N, we denote by Ck(ω) the set of all ∈ continuous real-valued functions u defined on ω such that ∂αu is continuous on ω for every multi-index α, α k. When k = 0, we shall write C(ω) in lieu of C0(ω). For k = , | | ≤ ∞ C∞(ω) will denote the intersection Ck(ω). k≥0 We shall also requirespaces of functions definedover the closure ω¯ of an open set ω Rd. T ⊂ For k N, Ck(ω¯) will signify the set of all u Ck(ω) such that ∂αu can be continuously ∈ ∈ extended from ω onto ω¯ for every multi-index α, α k. Further, we define C∞(ω¯) as the | | ≤ intersection Ck(ω¯). The notation C0(ω¯) is abbreviated to C(ω¯). k≥0 Assuming that ω is a bounded open set in Rn and k N, the linear space Ck(ω¯) is a T ∈ Banach space equipped with the norm u = max sup ∂αu(x) . k kCk(ω¯) |α|≤k x∈ω| | For k N we denote by Ck,1(ω¯) the set of all u Ck(ω¯) such that the quantity ∈ ∈ ∂αu(x) ∂αu(y) u = max sup | − | | |Ck,1(ω¯) |α|=k x6=y,x,y∈ω x y | − | is finite. Ck,1(ω¯) is a Banach space with the norm u = u + u . k kCk,1(ω¯) k kCk(ω¯) | |Ck,1(ω¯) Clearly, Ck+1(ω¯) Ck,1(ω¯). When u belongs to C0,1(ω¯), it is said to be Lipschitz continuous ⊂ on ω¯. The support, suppu, of a continuous function u defined on an open set ω is the closure in ω of the set x ω : u(x) = 0 ; in other words, suppu is the smallest closed subset of ω { ∈ 6 } such that u = 0 on ω suppu. For k = 0,1,..., , Ck(ω) denotes the set of all u Ck(ω) \ ∞ 0 ∈ whose support is a bounded (and, by definition, closed) subset of ω. 6 1.1.2 Spaces of integrable functions For p 1 and an open set ω Rn, L (ω) will denote the set of all real-valued Lebesgue p ≥ ⊂ measurable functions u defined on ω such that up is integrable on ω with respect to the | | Lebesgue measure dx = dx ... dx ; it is implicitly assumed that any two functions which 1 n are equal almost everywhere (i.e., equal, except maybe on a set of zero Lebesgue measure) are identified. L (ω) is a Banach space equipped with the norm p 1/p u = u(x)pdx . k kLp(ω) | | (cid:18)Zω (cid:19) When p = 2, L (ω) is a Hilbert space with the inner product 2 (u,v) = u(x)v(x) dx. ω Zω In the case when ω Ω, we write (, ) in lieu of (, ) . Ω ≡ · · · · L (ω) denotes the set of all real-valued Lebesgue measurable functions u defined on ω ∞ such that u has finite essential supremum; the essential supremum of u is defined as the | | | | infimum of the set of all positive real numbers M such that u M almost everywhere on | | ≤ ω. Again, any two functions that are equal almost everywhere on ω are identified. L (ω) is ∞ a Banach space with norm u = ess.sup u(x) . k kL∞(ω) x∈ω| | Ho¨lder’s Inequality. Let u L (ω) and v L (ω), where 1/p+1/q = 1, 1 p,q . Then p q ∈ ∈ ≤ ≤ ∞ uv L (ω) and 1 ∈ u(x)v(x) dx u v . ≤ k kLp(ω)k kLq(ω) (cid:12)Zω (cid:12) (cid:12) (cid:12) In the special case when p(cid:12)= q = 2, this (cid:12)inequality is referred to as the Cauchy-Schwarz (cid:12) (cid:12) Inequality. 1.1.3 Sobolev spaces Given that ω is an open set in Rn, k a non-negative integer and 1 p , we define the ≤ ≤ ∞ Sobolev space Wk(ω)= u L (ω) : ∂αu L (ω), α k , p { ∈ p ∈ p | | ≤ } and equip it with the Sobolev norm defined by 1/p ∂αu p , if 1 p < , |α|≤kk kLp(ω) ≤ ∞ u = k kWpk(ω)  (cid:16)P (cid:17)  u = max ∂αu , if p = . k kW∞k(ω) |α|≤kk kL∞(ω) ∞  The associated Sobolev seminorm is defined by 1/p ∂αu p , if 1 p < , |α|=kk kLp(ω) ≤ ∞ u = | |Wpk(ω)  (cid:16)P (cid:17)  max ∂αu , if p = . |α|=kk kL∞(ω) ∞   7 In these definitions the derivatives are to be understood in the sense of distributions. The Sobolev space Wk(ω) is a Banach space with the norm , 1 p , k 0. p k · kWpk(ω) ≤ ≤ ∞ ≥ Specifically, for p = 2, the normed linear space Wk(ω) is a Hilbert space with the inner 2 product (u,v) = (∂αu,∂αv) , Wk(ω) ω 2 |α|≤k X where (, ) denotes the inner product in L (ω). ω 2 · · Finersmoothnesspropertiesofintegrablefunctionscanbedetectedbyconsideringfractional- order Sobolev spaces. Given that s is a positive real number, s N, let us write s = m+σ, 6∈ where 0 < σ < 1 and m = [s] is the integer part of s. The fractional-order Sobolev space Ws(ω), 1 p < , is the set of all u Wm(ω) such that p ≤ ∞ ∈ p 1/p Dαu(x) Dαu(y)p |u|Wps(ω) =  | x −y n+σp | dxdy < ∞, |αX|=mZωZω | − |  with the usual modification when p = . The fractional-order Sobolev norm of index s is ∞ defined by 1/p p p u + u , if 1 p < , u = k kWpm(ω) | |Wps(ω) ≤ ∞ k kWps(ω)  n o   u + u , if p = . Wm(ω) Ws (ω) k k ∞ | | ∞ ∞  The fractional-order Sobolev space Ws(ω) is a Banach space with this norm. p When p =2 we shall write Hs in place of Ws to signify the fact that we are dealing with 2 a Hilbert space. We denote by Hs(ω) the closure of C∞(ω) in the norm of Hs(ω); when ω is 0 0 a Lipschitz domain and 1/2 <s < 3/2, this space coincides with the set of all those functions in Hs(ω) whose trace on ∂ω is equal to zero. 8 2 Motivation: Linear problems and adjoint equations In this section we present an overview of the general theoretical framework of adjoint–based a posteriori error estimation developed by C. Johnson and R. Rannacher and their collab- orators. For a detailed discussion, we refer to the series of articles [23, 36, 71, 79], and the references cited therein. To this end, we introduce the a posteriori error estimation for linear problemsinSection2.1. TheninSection2.2wegiveaframeworkforderivingadjointproblem for linear primal problems. This is then applied to the linear advection equation in Section 2.3. A numerical example in Section 2.4 highlights the practical importance of adjoint–based refinement. 2.1 Error estimation for linear problems We begin by considering a linear problem Lu = f in Ω, Bu= g on Γ, (1) where f L (Ω), g L (Γ), L denotes a linear differential operators on Ω, and B denotes a 2 2 ∈ ∈ linear boundary operator on Γ. Let the linear problem (1) be discretized as follows: Find u V such that h h,p ∈ (u ,v ) = ℓ(v ) v V , (2) h h h h h,p B ∀ ∈ whereV isafiniteelementspaceconsistingofpiecewisepolynomialfunctionsofdegreepon h,p a partition of the domain Ω in elements κ of size h. Furthermore, (, ) :V V R h h T ∈ T B · · × → is a bilinear form and l() : V R is a linear form including the forcing function f and the · → boundary value function g. Here, V is some suitably chosen function space including the analytical solution u V to the primal problem and satisfying V V. h,p ∈ ⊂ Furthermore, let us assume that the discretization (2) is consistent, i.e., the analytical solution u V satisfies the following equation: ∈ (u,v) = ℓ(v) v V. (3) B ∀ ∈ Inmanyproblemsofphysicalimportancethequantitiesofinterestmaybeaseriesoftarget or error functionals J (), i = 1,...,N, N 1, of the solution. Relevant examples include i · ≥ the mean value of the solution, the mean flow across a line, the point value of the solution or different scalar quantities which can be computed from the solution u. For compressible flows, which are not covered by the theory in this introductory section, such a quantity J(u) could represent an aerodynamic force coefficient, like the drag, lift or moment coefficient. For simplicity, werestrictourselvestothecaseofasinglelineartargetfunctional, i.e., N = 1, and write J() J (); for the extension of the proceeding theory to multiple target functionals, 1 · ≡ · see Section 4.3; cf., also, [60]. In order to obtain a computable a posteriori bound on the error between the true value of the functional J(u) and the computed value J(u ), we begin h by noting the Galerkin orthogonality of the discretization (2): (u,v ) (u ,v )= (u u ,v )= 0 v V . (4) h h h h h h h,p B −B B − ∀ ∈ This will be a key ingredient in the following a posteriori error analysis. 9 We now introduce the following adjoint problem: find z V such that ∈ (w,z) = J(w) w V; (5) B ∀ ∈ We assume that (5) possesses a unique solution; clearly, the validity of this assumption depends on both the definition of (, ) and the choice of the functional under consideration. B · · Important examples which are covered by our hypothesis are discussed below, cf. [77]. For the proceeding error analysis, we must therefore assume that the adjoint problem (5) is well–posed. Under this assumption, employing the Galerkin orthogonality property (4) we deduce the following error representation formula: J(u) J(u ) = J(u u ) = (u u ,z) h h h − − B − = (u u ,z z ) h h B − − = ℓ(z z ) (u ,z z ) (6) h h h − −B − for all z in the finite element space V . On the basis of the general error representation h h,p formula (6), a posteriori estimates which provide upper bounds on the true error in the computedtargetfunctionalJ()maybededuced. Thesimplestapproachistofirstdecompose · the right–hand side of (6) as a summation of local error indicators η over the elements κ in κ the computational mesh , i.e., we write h T J(u) J(u ) = ℓ(z z ) (u ,z z ) (u ,z z ) = η ; (7) h h h h h h κ − − −B − ≡ R − κX∈Th then,uponapplicationofthetriangleinequality,wededucethefollowingweightedaposteriori error bound. Theorem 2.1 Let u and u denote the solutions of (1) and (2), respectively, and suppose h that the adjoint problem (5) is well–posed. Then, the following aposteriori error bound holds: J(u) J(u ) (u ,z z ) η (8) h |Ω| h h κ | − | ≤ R − ≡ | | κX∈Th for all z in V . h h,p We remark that the local error indicators η appearing on the right–hand side of (9) in- κ volve the multiplication of finite element residuals depending only on u with local weighting h terms involving the difference between the adjoint solution z satisfying (5) and its projec- tion/interpolant z onto the finiteelement space V . These weights representthe sensitivity h h,p of theerrorin thetarget functionalJ() withrespecttovariations of thelocal element residu- · als; indeed, they provide invaluable information concerning the global transport of the error, which is essential for efficient error control. Since the solution to the adjoint problem is usually unknown analytically it must be numerically approximated, cf. [23, 43, 58]. Replacing the unknown exact adjoint solution z in (7) by a numerical approximation z¯ / V , we obtain following approximate error h h,p ∈ representation J(u) J(u )= (u ,z z ) (u ,z¯ z ) = η¯ . (9) h h h h h h κ − R − ≈ R − κX∈Th Note that the so-called adjoint-based indicators η in (9) can be used to drive an adaptive κ algorithm targeted at the accurate and efficient approximation of the target quantity J(u). In the following sections we give some examples of adjoint problems. 10