- Email: [email protected]
The effect of stabilization in finite element methods for the optimal boundary control of the Oseen equations
Finite Elements in Analysis and Design 41 (2004) 229 – 251 www.elsevier.com/locate/finel
The effect of stabilization in finite element methods for the optimal boundary control of the Oseen equations Feby Abrahama , Marek Behrb,∗ , Matthias Heinkenschlossc a Department of Mechanical Engineering and Materials Science, Rice University MS 321, 6100 Main Street, Houston,
TX 77005, USA b Lehrstuhl für Numerische Mechanik, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany c Department of Computational and Applied Mathematics, Rice University MS 134, 6100 Main Street, Houston, TX 77005, USA
Received 26 May 2003; received in revised form 14 May 2004; accepted 15 June 2004
Abstract We study the effect of the Galerkin/Least-Squares (GLS) stabilization on the finite element discretization of optimal control problems governed by the linear Oseen equations. Control is applied in the form of suction or blowing on part of the boundary. Two ways of including the GLS stabilization into the discretization of the optimal control problem are discussed. In one case the optimal control problem is first discretized and the resulting finitedimensional problem is then solved. In the other case, the optimality conditions are first formulated on the differential equation level and are then discretized. Both approaches lead to different discrete adjoint equations and, depending on the choice of the stabilization parameters and grid size, may significantly affect the computed control. The effect of the order in which the discretization is applied and the choice of the stabilization parameters are illustrated using two test problems. The cause of the differences in the computed controls are explored numerically. Diagnostics are introduced that may guide the selection of sensible stabilization parameters. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Optimal boundary control; Stabilized finite element methods; Oseen equations; Solution accuracy
1. Introduction Stabilized finite element methods (FEMs) are frequently and successfully used to discretize advection-dominated partial differential equations (PDEs) [1], or to circumvent the compatibility ∗ Corresponding author.
E-mail addresses: [email protected] (F. Abraham), [email protected] (M. Behr), [email protected] (M. Heinkenschloss). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.06.001
230
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
conditions restricting the choice of interpolation function spaces [2]. One question that arises in the application of these methods is the choice of the stabilization parameter. For many PDE model problems this issue has been studied analytically and numerically [3,4]. However, the impact on the choice of the stabilization parameter when the stabilized FEM is used in the context of optimal control is not well studied. It is not true that a scheme which gives good approximations to the PDE solution for a fixed simulation is also guaranteed to provide good approximations in the context of optimal control problems. Solutions of optimal control problems are characterized by the original governing PDE as well as another PDE—the so-called adjoint equation. Depending on the approach chosen, a discretization of the original governing PDE may imply a discretization scheme for the adjoint equation with poor approximation properties. The present paper investigates this issue for a class of boundary control problems discretized using the Galerkin/least-squares (GLS) approach [5]. The purpose of this paper is to numerically investigate the influence of stabilization parameters on the solution of optimal control problems related to the optimal boundary control of the Navier–Stokes equations. We will demonstrate that the way the stabilization is introduced into the formulation of the optimal control problem, and the choice of stabilization parameters, can have a significant impact on the computed control. Furthermore, we point out some diagnostic tools that may help to assess the quality of the computed control and guide the choice of stabilization parameters. We consider a class of linear quadratic optimal control problems governed by the Oseen equations, and we study the effect of the GLS-stabilized finite element method on the computed control. The linear Oseen equations were chosen as the governing state equation instead of the nonlinear Navier–Stokes equations, because the resulting optimal control problem has a unique solution and the first-order optimality conditions are necessary and sufficient. Optimal control problems governed by the nonlinear Navier–Stokes equations may have local solutions and their solution requires iterative methods. Since it is difficult to separate the possible effects of local solutions and iterative solvers on the computed optimal control from the effects of the stabilization on the computed optimal control, we have chosen the Oseen equations. However, the linear quadratic optimal control problems studied in this paper are closely related to the subproblems that arise in the solution of boundary control problems governed by the Navier–Stokes equations using Newton or sequential quadratic programming methods (see, e.g., [6,7]). Consequently, the results reported in this paper are also relevant for the optimal control of Navier–Stokes flow. A related paper [8] studies the effect of the streamline-upwind/Petrov–Galerkin stabilized finite element method on the numerical solution of linear quadratic distributed optimal control problems governed by an advection-diffusion equation. That paper contains both analytical results that describe the convergence of the computed control (state/adjoint) to the exact control (state/adjoint) as the grid is refined, as well as numerical convergence studies for the simpler optimal control problem. The analytical results in [8] are accurate asymptotically, but do not describe the numerical behavior well for first-order finite elements on relatively coarse, but for practical purposes acceptable grids. Therefore, the present paper focuses on a numerical study. We expect that most of the analytical results in [8] can be extended to the problems and the GLS stabilization method considered in this paper. However, such a theoretical treatment is beyond the scope of this paper. We consider a viscous incompressible fluid occupying a bounded region ⊂ Rnsd , where nsd is the number of space dimensions. The symbols u and p represent the velocity and pressure. The boundary * of is decomposed into three disjoint segments h , d and g . Suction and blowing control is applied
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
on g . We are interested in the solution of the following problem: minimize 2 J (u, g) = J1 (u) + |g|2 dx 2 g
231
(1)
subject to (a · ∇)u − ∇ · (u, p) = f ∇ ·u=0
in ,
(2a)
in ,
(u, p) · n = h
(2b) on h ,
(2c)
u=d
on d ,
(2d)
u=g
on g .
(2e)
The above equations represent the momentum and continuity equations subject to Neumann- and Dirichlettype boundary conditions. In (2), the functions f, h and d are given, and the control g has to be determined as the solution of the optimization problem. The regularization parameter 2 > 0 is given. The stress tensor is given by (u, p) = −pI + 2(u),
(u) =
1 (∇u + ∇uT ), 2
(3)
where is dynamic viscosity and I denotes the identity tensor. We consider two objective functions J (u, g) with 1 J1 (u) = |∇ × u|2 dx (4) 2 ∗ and
J1 (u) = 2
∗
(u): (u) dx,
(5)
where ∗ ⊂ . These objective functions—vorticity magnitude and dissipation due to viscous stress—represent design goals encountered in biomedical engineering. For example, artificial flow devices such as blood pumps induce non-physiological levels of shear stress, leading to blood damage [9] or thrombus formation (clotting) [10]. We give a precise statement of the optimal control problem in Section 2.1, as well as establish existence of a unique solution and formulate the first-order necessary and sufficient optimality conditions. In Section 2.2 we recall the GLS-stabilized finite element method for the discretization of the state equation. The need for stabilization in (2) is twofold. Firstly, GLS removes the under-diffusivity of Galerkin methods which leads to oscillatory solutions in advection-dominated flows. Secondly, GLS circumvents the Babuška–Brezzi (inf–sup) condition, and allows the use of equal-order finite elements for velocities and pressure. For the numerical solution of linear quadratic optimal control problems, there are at least two approaches. In the first approach, called discretize-then-optimize (do), the objective function and governing
232
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
PDEs are discretized. In our case, the Oseen equation is discretized using the GLS-stabilized finite element formulation. This leads to a large-scale finite-dimensional quadratic programming problem, or equivalently, to a large-scale linear system that is solved using suitable numerical linear algebra tools. In the second approach, called optimize-then-discretize (od), the first-order necessary and sufficient optimality conditions are formed and then discretized. The first-order optimality conditions consist of the Oseen equations, the so-called adjoint PDEs—which have a structure similar to the Oseen equations—and an algebraic equation that links controls and adjoint variables. These PDEs are then individually discretized, in our case using the GLS formulation applied to the Oseen and the adjoint PDEs. This also results in a large-scale linear system. These two approaches will be described in Sections 2.3 and 2.4, respectively. The linear systems arising in either approach are not the same, due to the way the stabilization terms affect the adjoint PDEs. These differences and the choice of the GLS stabilization parameter will be explored numerically in Section 4. We will see that an inappropriate choice of the stabilization parameter leads to significant differences in the computed controls. Moreover, the computed control can be very sensitive to a scalar weighting parameter in the stabilization term. In our results, the solutions computed by the optimize-then-discretize approach are more sensitive to the choice of the stabilization parameter and the scalar weighting term. If the stabilization is computed using an element length based on the direction of the advective field a, the optimal controls computed by either approach are in good agreement on fine grids. 2. A model problem 2.1. Formulation of the optimal control problem In the following, we use the spaces L2 () and H 1 (), which are defined in the usual way [11], and V = {v ∈ H 1 () | v = 0 on d }. We set H1 () = [H 1 ()]nsd , V = V nsd , L2 (h ) = [L2 (h )]nsd , and L2 (g ) = [L2 (g )]nsd . The Dirichlet boundary conditions (2d) and (2e) can be implemented through interpolation [12], weakly through a Lagrange multiplier technique [13], or via a penalty approach [7,14]. We use interpolation to implement the Dirichlet boundary conditions (2d) with fixed data, and replace (2e) by the penalized Neumann boundary condition (u, p) · n +
1
u=
1
g
on g ,
(6)
where > 0 is a given penalty parameter. This choice enables us to look for controls g in L2 (g ), instead of H1/2 (g ). The weak form of the partial differential equation (2a)–(2d), (6) is given as follows: Find u ∈ H1 () with u = d on d and p ∈ L2 (), such that 1 (a · ∇)u · v dx + (v): (u, p) dx − q∇ · u dx + v · (u − g) dx g = v · h dx + f · v dx (7) h
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
233
for all v ∈ V and all q ∈ L2 (). It is also possible to implement (2d) via a penalized Neumann approach. See Remark 2.3 below. Eq. (7) motivates the definition of the bilinear forms a(u, v) = 2 (v): (u) dx ∀v, u ∈ H1 (), (8a)
c(u, v) =
(a · ∇)u · v dx
∀v, u ∈ H1 (),
(8b)
b(u, q) =
q∇ · u dx
∀u ∈ H1 (), ∀q ∈ L2 (),
(8c)
h · v dx
∀h ∈ L2 (h ), ∀v ∈ V,
(8d)
g · v dx
∀g ∈ L2 (h ), ∀v ∈ V
(8e)
h, vh = g, vg =
h
g
and of the linear functional f, v = f · v dx ∀v ∈ V.
(8f)
For a given > 0, we consider the optimal control problem 2 minimize J1 (u) + |g|2 dx, 2 g subject to
(9a)
a(u, v) + c(u, v) − b(v, p) − b(u, q) + −1 u, vg − −1 g, vg = h, vh + f, v
∀v ∈ V, ∀q ∈ L2 ().
(9b)
Lemma 2.1. Let > 0 be given and assume that there exists ud ∈ H1 () with ud =d on d . If a ∈ H1 (h ) satisfies ∇ · a = 0 and a · n 0 on \d , then for any h ∈ L2 (h ), g ∈ L2 (g ) and f ∈ L2 () Eq. (9b) has a unique solution (u, p) ∈ H1 () × L2 (). Moreover, there exists a constant C > 0 (dependent on , but independent of h, g, f and ud ) such that uH1 () + pL2 () C(hL2 (h ) + gL2 (h ) + fL2 () + ud H1 () ). Proof. The bilinear form b satisfies the inf–sup condition inf
sup b(v, q)/(qL2 () vV ) > 0
q∈L2 () v∈V
[15, p. 81]. Integration by parts yields 1 1 2 c(v, v) = − ∇ · a |v| dx + a · n |v|2 dx 2 2 * With the assumptions on a, this implies c(v, v) 0
∀v ∈ V.
∀v ∈ H1 ().
(10)
234
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Hence, there exists > 0 such that a(v, v) + c(v, v) + −1 v, vg v2H1 ()
∀v ∈ V.
(11)
Continuity of the linear and bilinear forms (8) follows using the arguments in [14]. The assertion now follows from the theory of saddle point problems [16, Section II.1], [15, Section 4]. Theorem 2.2. Let the assumptions of Lemma 2.1 be satisfied. If J1 is weakly lower semicontinuous on V, then (9) admits a solution (u , p , g ) ∈ H1 () × L2 () × L2 (g ). If J1 is convex, e.g., if J1 is given by (4) or (5), the solution is unique. Proof. The result follows from standard arguments, see, e.g., the proof of Theorem 3.5 in [14]. Remark 2.3. It is possible to also implement the Dirichlet boundary condition (2d) via a penalized Neumann condition 1 1 1 (u, p) · n − (a · n) u + u = d on d . (12) 2 The weak form of the partial differential equation (2a)–(2c), (6), (12) is given as follows: Find u ∈ H1 () and p ∈ L2 (), such that 1 (a · ∇)u · v dx + (v): (u, p) dx − q∇ · u dx − (a · n)|u|2 dx 2 d 1 1 v · (u − g) dx + v · (u − d) dx + g d = v · h dx + f · v dx (13) h
for all v ∈ and all q ∈ L2 (). Lemma 2.1 with V replaced by H1 () still holds, since (10) and the inclusion of the term − 21 (a · n) u in (12) imply H1 ()
a(v, v) + c(v, v) + −1 v, vg v2H1 ()
∀v ∈ H1 ().
(14)
The convergence behavior of the solutions (u , p , g ) of the optimal control problem (9) with penalized Neumann boundary control to the solution of (1) and (2) as → 0 is discussed in [14] in the context of Navier–Stokes equations. For the purpose of this study, we view (9) with a small but fixed > 0 as our optimal control problem. The Lagrangian associated with the optimal control problem (9) is given by 2 L(u, p, g, , )=J1 (u) + |g|2 dx 2 g + a(u, ) + c(u, ) − b(, p) − b(u, ) + −1 u, g − −1 g, g − h, h − f, . (15) Lemma 2.1 ensures that (9) satisfies a constraint qualification. Hence, the necessary and sufficient optimality conditions for (9) are obtained by setting the Fréchet derivatives of the Lagrangian with
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
235
respect to the adjoint variables (, ), with respect to the state variables (u, p) and with respect to the controls g equal to zero. The necessary and sufficient optimality conditions consist of the state equation a(u, v) + c(u, v) − b(v, p) − b(u, q) + −1 u, vg − −1 g, vg = h, vh + f, v
∀v ∈ V, ∀q ∈ L2 (),
(16a)
the adjoint equation a(v, ) + c(v, ) − b(, q) − b(v, ) + −1 v, g = −Du J1 (u), v
∀v ∈ V, ∀q ∈ L2 ()
(16b)
∀z ∈ L2 (g ).
(16c)
and the gradient equation 2 g, zg =
1
, zg
In the adjoint equation (16b), Du J1 (u) denotes the Fréchet derivative of J1 with respect to u and ·, · is the duality pairing between (H1 ())∗ —the dual of H1 ()—and H1 (), i.e., Du J1 (u), v is the application of the Fréchet derivative of J1 to the function v. In the case of our objective functions (4) and (5) these are simply given as Du J1 (u), v = (∇ × u) · (∇ × v) dx (17) ∗
and
Du J1 (u), v = 4
∗
(u): (v) dx,
(18)
respectively. The adjoint equation (16b) may formally be interpreted as the weak form of the adjoint partial differential equation −(a · ∇) − (∇ · a) − ∇ · (, ) = −Du J1 (u) ∇ ·=0
in ,
1 2
on g .
(19c) (19d)
1
and Eq. (16c) simply states that g=
on h ,
on d ,
(, ) · n + (a · n) +
(19a) (19b)
(, ) · n + (a · n) = 0 =0
in ,
=0
on g ,
(19e)
236
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
2.2. Discretization of the state equations For the discretization of the state equation we apply the GLS-stabilized finite element method using piece-wise-linear polynomials for both the velocity and the pressure. We divide our domain into nel subdomains e , e = 1, 2, . . . , nel . We assume that our triangulation is such that the controlled boundary g can be written as the union of edges or faces of elements e , e = 1, 2, . . . , nel . We use piecewiselinear polynomials to approximate both the velocities and pressures. Hence, the discretized velocities are in H1h () = [H 1h ()]nsd , where H 1h () = { h ∈ C 0 () | h |e ∈ P 1 , e = 1, 2, . . . , nel } and the discretized pressures are in L2h () = { h ∈ C 0 () | h |e ∈ P 1 , e = 1, 2, . . . , nel }. Furthermore, we define V h = { h ∈ H 1h () | h = 0 on d }, Vh = [V h ]nsd , L2h (g ) = { h |g | h ∈ H 1h ()} and L2h (g ) = [L2h (g )]nsd . We let Ih d denote the piecewise linear interpolant of the Dirichlet boundary data d. The GLS-stabilized finite element formulation for the state equation (2a)–(2c), (6) is given as follows: Find uh ∈ H1h () with uh = Ih d on d and ph ∈ L2h () such that h h v · (a · ∇)u dx + (vh ) : (uh , ph ) dx
+
nel e
e [(a · ∇)vh − ∇ · (vh , q h )] · [(a · ∇)uh − ∇ · (uh , p h )] dx
e=1 1 h h + q ∇ · u dx + =
h
vh · h dx +
g
vh · (uh − gh ) dx
f · vh dx +
nel e=1
e
e [(a · ∇)vh − ∇ · (vh , q h )] · f dx
(20)
for all vh ∈ Vh and all q h ∈ L2h (). Our choice [3] for the stabilization parameter e is
e =
2|a| he
2
+
4 h2e
2 −1/2 ,
(21)
where is a factor typically taken to be 1.0 and the element length he is computed using one of two choices diag and adv . Let h∗e be the longest edge in an element. The diag element length is given by he = h∗e . The adv element length, introduced in [17], is defined as ∗ if |aeh | = 0, he
−1 he = (22) otherwise, a |se · ∇Na | where aeh is the interpolated advection velocity at the center of element e, se = aeh /|aeh |, and Na represents the basis functions. Both these element length definitions are commonly used in the computation of
e . Notice that the adv element length is identical to the diag element length in regions where there is stagnation in the advection field. In many flow scenarios with extended and frequent regions of stagnation, the diag element length is used to avoid abrupt switching of the element length type from streamline interpolation to longest edge of the element, as is the case in the adv element length definition.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
237
2.3. Discretization of the optimization problem A frequently used approach for the numerical solution of an optimal control problem is to discretize the optimal control problem first and then solve the resulting finite-dimensional nonlinear (in our case quadratic) programming problem using a suitable optimization algorithm. We refer to this approach as the discretize-then-optimize (do) approach. Using the GLS-stabilized finite element method described in Section 2.2 for the discretization of the state equation, our discretization of the optimal control problem (9) is given by 2 J1 (uh ) + |gh |2 dx, 2 g (20), uh ∈ H1h (), ph ∈ L2h (), gh ∈ L2h (g ), uh = Ih d on d .
minimize subject to with
(23)
The Lagrangian for the discretized problem (23) is given by Lh (uh , ph , gh , h , h ) 2 h = J1 (u ) + |gh |2 dx 2 g h h + · (a · ∇)u dx + (h ): (uh , ph ) dx
+
e=1 + −
nel
e
h
e [(a · ∇)h − ∇ · (h , h )] · [(a · ∇)uh − ∇ · (uh , p h )] dx h
∇ · u dx +
h
h · h dx −
1
g
h · (uh − gh ) dx
f · h dx −
nel e=1
e
e [(a · ∇)h − ∇ · (h , h )] · f dx.
(24)
The necessary and sufficient optimality conditions for (23) are obtained by setting the derivatives of the Lagrangian with respect to the discrete adjoint variables (h , h ), the discretized state variables (uh , ph ) and the discretized controls gh to zero. Setting the derivative of Lh with respect to the discrete adjoint variables (h , h ) to zero gives the discretized state equation (20). Setting the derivative of Lh with respect to the discretized state variables (uh , ph ) to zero gives the discrete adjoint equation h h · (a · ∇)v dx + (h ): (vh , q h ) dx
+
e=1 +
nel
e
h
e [(a · ∇)h − ∇ · (h , h )] · [(a · ∇)vh − ∇ · (vh , q h )] dx h
∇ · v dx +
= −Du J1 (uh ), vh
1
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(25)
238
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
where Du J1 (uh ), vh is defined as in (17) or (18), respectively. Finally, setting the derivative of Lh with respect to the discretized control variables gh to zero gives the discrete gradient equation 1 2 h h g · z dx = h · zh dx ∀zh ∈ L2h (g ). (26)
g
g
The necessary and sufficient conditions (25), (26) and (20) lead to a system of linear equations
Hdo 0 A
0 GT −B
AT 0 u do T gdo = 0 . −B f 0 do
(27)
Here u do and do are the vectors containing the values of (uh , ph ) and (h , h ) at the grid points, respectively. 2.4. Discretization of the optimality conditions The optimality conditions (16) are necessary and sufficient. Therefore, we can tackle the optimality conditions directly to compute an approximate solution of our optimal control problem. This is called the optimize-then-discretize (od) approach. This approach requires us to discretize the state and the adjoint PDE (19). In principle, these discretizations do not have to be the same. However, we will use the same mesh and the GLS-stabilized finite element method for the discretization of both. We use the same notation as in Section 2.2. Our discretization of the state equation is given in (20). The GLS-stabilized finite element method applied to the adjoint equation (19), yields h h · (a · ∇)v dx + (vh ): (h , h ) dx
+
nel e
e [−(a · ∇)vh − (∇ · a)vh − ∇ · (vh , q h )] · [−(a · ∇)h − (∇ · a)h − ∇ · (h , h )] dx
e=1 1 h h + q ∇ · dx +
= −Du J1 (uh ), vh s
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(28)
where Du J1 (uh ), vh s is given as follows. For the objective function (4), Du J1 (uh ), vh s = (∇ × uh ) · (∇ × vh ) dx ∗ nel
+
e=1
e ∩∗
e (∇ × uh ) · (∇ × [−(a · ∇)vh
− (∇ · a)vh − ∇ · (vh , q h )]) dx,
(29)
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
239
and for the objective function (5), Du J1 (uh ), vh s =4 (uh ): (vh ) dx ∗
+ 4
nel e=1
e ∩∗
e (uh ): (−(a · ∇)vh
− (∇ · a)vh − ∇ · (vh , q h )) dx.
(30)
The discretization of (16c) is again given by (26). The necessary and sufficient conditions (20), (28) and (26) lead to a system of linear equations 0 A 0 Hod u od T T 0 (31) G −B god = 0 . f A −B 0 od Here u od and od are the vectors containing the values of (uh , ph ) and (h , h ) at the grid points, respectively. Note that because of the stabilization terms in (29) and in (30), the submatrix Hdo in (27) is different from the submatrix Hod in (31). Moreover, because the advective terms in adjoint equation (19) are different from those in the Oseen equation (2), the discrete adjoint equations (25) are different from the discretized adjoint equations (28). Hence the submatrix AT in (27) is also different from submatrix A in (31). 3. Implementation 3.1. Stabilization terms Our choice of piecewise-linear functions makes our discrete approximation a low-order method, i.e. the order of the function space is lower than the order of the highest derivative in (2), the stabilized finite element equivalent of which is (20). Therefore, the terms ∇ · (vh , q h ) and ∇ · (uh , ph ) in (20) are given by ∇q h and ∇p h , respectively. Dropping the term ∇ · 2(uh ) guarantees us only a weak consistency [18], since this term itself vanishes with refinement, as e approaches zero based on the relation in (21). Normally, a variational reconstruction of this term is done to achieve a “stronger consistency” for such lower-order choice of function space [18]. However, we have not done so here, to avoid the potential nonlinearity thus introduced. A similar inconsistency is introduced into (25) and (28). The implemented version of (25) is given by h h · (a · ∇)v dx + (h ): (vh , q h ) dx
+
e=1 +
nel
e
h
e [(a · ∇)vh + ∇q h ] · [(a · ∇)h + ∇ h ] dx h
∇ · v dx +
= −Du J1 (uh ), vh
1
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(32)
240
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
and the implemented version of (28) is given by h h · (a · ∇)v dx + (vh ): (h , h ) dx
+
nel e
e [−(a · ∇)vh − (∇ · a)vh + ∇q h ] · [−(a · ∇)h − (∇ · a)h + ∇ h ] dx
e=1 1 h h + q ∇ · dx +
= −Du J1 (uh ), vh s where h
h
Du J1 (u ), v s =
∗
+
h
h
(33)
(∇ × uh ) · (∇ × vh ) dx
nel
Du J1 (u ), v s =4
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
e=1
and
g h
∗
+ 4
e ∩∗
e (∇ × uh ) · (∇ × [−(a · ∇)vh − (∇ · a)vh ]) dx,
(34)
(uh ): (vh ) dx
nel e=1
e ∩∗
e (uh ): (−(a · ∇)vh − (∇ · a)vh ) dx,
(35)
for the objective functions (4) and (5), respectively. 3.2. Differences The use of the GLS-stabilized finite element method in the discretization of the optimal control problem creates differences between the discretize-then-optimize approach and the optimize-then-discretize approach. Specifically, there are differences between the discrete adjoint equation (32) and the discretized adjoint equation (33). To explore how these differences impact the computed solution, we also implement variations of (33). 1. The left- and right-hand side in (33) contain ∇ · a terms, while the divergence of the advective field does not enter into the discrete adjoint equation (32). Therefore, we also compute the optimal controls using (20), (26) and (33), where the ∇ · a terms in (33) are dropped. The optimal controls computed this way will be labeled as od1. 2. If we compare (32) and (33) with ∇ · a terms in (33) replaced by zero, we see that the right-hand sides of the discretized adjoint equations (33) contain −(a · ∇)vh terms that are not present in the discrete adjoint equations (32). Therefore, we also compute the optimal controls using (20), (26) and (33) where all ∇ · a terms in (33) as well as the −(a · ∇)vh term on the right-hand side of (33) are dropped. The optimal controls computed this way will be labeled as od2.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
241
3. After the modifications of the discretized adjoint (33) described above in 2. have been performed, the only remaining difference between (32) and (33) is the sign of (a · ∇)vh and (a · ∇)h .
3.3. Sensitivity analysis We are interested in the sensitivity of the computed solution with respect to the stabilization factor . Eqs. (25), (26) and (20) reveal that the submatrix A and the right hand side component f in (27) depend on . Hence we write A( ) and f( ) instead of A and f, respectively. Moreover, we denote the solution of (27) by u do ( ), gdo ( ), do ( ). Examination of (25), (26) and (20) reveals that the submatrices A( ) and Hdo ( ) in (27) are of the form A( ) = A1 + A2 and Hdo ( ) = H1 , where we use subscript 1 for all those terms that do not depend on stabilization, and subscript 2 for those terms that are dependent on ( )=0. stabilization. Hence, if we use to denote differentiation with respect to , then A ( )=A2 and Hdo By the implicit function theorem, the derivatives of the solution of(27) satisfy
Hdo
0
A( )T
0 GT
A( ) −B
( ) u do
−BT gdo ( ) 0 do ( )
0
= − 0 A ( )
0
(A ( ))T
0
0
0
0
0 gdo ( ) + 0 , f ( ) do ( ) u do ( )
where u do ( ), gdo ( ) and do ( ) solve (27). The sensitivity of the solution u od ( ), god ( ) and od ( ) of (31) can be computed analogously. Note that in this case Hod also depends on and is of the form Hod ( ) = H1 + H2 , with Hod ( ) = H2 .
4. Numerical results In this section we study the effect of the GLS stabilization on the computed control for two test cases derived from commonly used model problems. In both cases, the advective field a is computed by solving the Navier–Stokes equations (a · ∇)a − ∇ · (a, p) = 0
∇ · a=0
on ,
on , (36)
with appropriate boundary conditions. The Navier–Stokes equations (36) are discretized using the GLSstabilized finite element method and the resulting nonlinear system is solved using Newton’s method with a residual tolerance of 10−16 . Note that while the “exact” advective field a is divergence-free, this is not true for the computed advective field, which is used as the coefficient in our computations. We solve (9) with parameters = 10−5 and 2 = 10−5 . As stated before, we use the notation do to refer to the control computed using the discretize-then-optimize approach (26), (20) and (32), and od, to refer to the control computed by the optimize-then-discretize approach (26), (20) and (33), respectively.
242
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Γt Γin Γg
Ω
Ω*
Γo
Γb
Fig. 1. Computational domain for test case 1.
Fig. 2. Finite element mesh for test case 1, discretization 1.
4.1. Test case 1 The first problem is modeled after the backward-facing step problem [7, pp. 1767,1769]. A schematic of the geometry is given in Fig. 1. The height of the inflow boundary is 0.5 and that of the outflow boundary is 1. The length of the narrower section of the channel is 1 and that of the wider section is 7, with the total horizontal length being 8. The advective velocity field a is computed by solving the Navier–Stokes equations (36) with the following boundary conditions. The inflow velocity is assumed parabolic with a profile a(0, y) = (8(y − 0.5)(1 − y), 0)T . At the outflow boundary o we impose traction-free boundary conditions for a in the x-direction. No-slip conditions are imposed at all other boundaries. We define the Reynolds number Re = Umax H /, where Umax is the maximum inlet velocity and H is the channel height. The Reynolds number for computing a is Re = 200. The boundary conditions for the Oseen equation are as follows. The inflow velocity field, different from that of a, is taken as parabolic given by u(0, y) = ((y − 0.5)(1 − y), 0)T . A traction-free boundary condition for the x-velocity is imposed at o . The controlled boundary g is the backward facing step. We allow suction and blowing in the normal direction. At the remaining boundaries b and t no-slip conditions are imposed for u. The objective function is given by (4), where ∗ is the subdomain of length 2 and height 0.5, as indicated in Fig. 1. We present results for two dynamic viscosity coefficients = 0.0005 and 0.00005, and for three discretizations. Discretization 1 consists of 1800 triangular elements as shown in Fig. 2. Discretization 2 consisting of 7200 triangular elements and is a uniform refinement of discretization 1. Discretization 3 is an uniform refinement of discretization 2 and consists of 28,800 elements. Fig. 3 shows the control velocities obtained using the various approaches considered here. We observe that the od and do approaches with the adv element length definition converge with mesh refinement. When we use the diag element length, we observe that the discretize-then-optimize approach converges to the solution obtained using the adv element length. However, we can see that there is a strong oscillatory behavior in the controls for optimize-then-discretize diag case, especially for the smallest viscosity, i.e. = 0.00005. Fig. 4 shows the dependence of control velocities on the factor , with increasing oscillatory behavior as increases. Fig. 5 shows the sensitivity of the controls to . The sensitivity of the computed control with respect to is large when the diag element length definition is chosen. In particular for the od approach with diag element length definition large sensitivities are observed near the corner of
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Controlling velocity for µ = 0.0005
0.5
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
Controlling velocity for µ = 0.00005
0.5
0.45
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
do od do od
0.1 0.05
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015 0.02 X velocity
0.025
0.03
0.035
0
0.005 0.01
0.5
0.45
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
0.5
0.3 0.25 0.2 0.15 do od do od
0.05
0.035 0.04 0.045
0.3 0.25 0.2 0.15
0.1
0.015 0.02 0.025 0.03 X velocity
Controlling velocity for µ = 0.00005
Controlling velocity for µ = 0.0005
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015 0.02 X velocity
0.025
0.03
0.035
0
Controlling velocity for µ = 0.0005
0.5
0.005
0.01
0.015
0.02 0.025 X velocity
0.03
0.035
0.04
0.045
0.04
0.045
Controlling velocity for µ = 0.00005
0.5
0.45
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
243
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
do od do od
0.1 0.05
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015
0.02
X velocity
0.025
0.03
0.035
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
X velocity
Fig. 3. The control velocities obtained for the two optimize-then-discretize and discretize-then-optimize approaches using the diag and adv element length definition with = 1 for test case 1. The ordering of these plots is = 0.0005 on the left and = 0.00005 on the right with discretizations 1–3 from top to bottom.
the backward facing step. This correlates with the ‘spikes’ in the computed control for this case shown in Fig. 3. The size of the sensitivities are relatively small for both, the discretize-then-optimize and the optimize-then-discretize approach when the adv element length definition is chosen. Again this correlates with the observed grid convergence for both approaches in case of the adv element length definition (see Fig. 3).
244
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Controlling velocity for µ = 0.00005
0.5
0.4
0.35
0.35 Y co-ordinate
Y co-ordinate
0.45
0.4
0.3 0.25 0.2 0.15 0.1 0.05
Controlling velocity for µ = 0.00005
0.5
0.45
od od od od
0.3 0.25 0.2 0.15
diag (α = 0.5) diag (α = 1.0) diag (α = 2.0) adv (α = 1.0)
do do do do
0.1 0.05
0 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 X velocity
diag (α = 0.5) diag (α = 1.0) diag (α = 2.0) adv (α = 1.0)
0 0.1
0
0.005
0.01
0.015
0.02 0.025 X velocity
0.03
0.035
0.04
Fig. 4. The control velocities obtained using different choices of using diag element length definition for = 0.00005 in discretization 3 for test case 1. The plot on the left uses optimize-then-discretize approach and the one on right uses discretize-then-optimize approach. For comparison we also present adv element length definition with = 1.
Sensitivity plot for µ = 0.00005, α = 0.5
0.5 0.45
0.45
0.4
0.4
0.35
0.35 Y co-ordinate
Y co-ordinate
Sensitivity plot for µ = 0.00005, α = 1.0
0.5
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
0.1
0.1
do diag od diag do adv od adv
0.05
do diag od diag do adv od adv
0.05
0
0 0
0.1
0.2
0.3 0.4 0.5 Weighted sensitivity
0
0.7
0.6
0.2
0.4
0.6 0.8 Weighted sensitivity
1
1.2
1.4
Sensitivity plot for µ = 0.00005, α = 2.0
0.5 0.45 0.4
Y co-ordinate
0.35 0.3 0.25 0.2 0.15 0.1
do diag od diag do adv od adv
0.05 0 0
1
2
3
4
5
6
Weighted sensitivity
Fig. 5. The weighted sensitivity |(g ( ))i |/|g( )i | of the controls for the discretize-then-optimize and optimize-then-discretize approaches using the diag and adv element length definitions for = 0.00005 and discretization 3 for test case 1, with = 0.5, = 1.0 and = 2.0.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
245
Controlling velocity for µ = 0.00005 0.5 0.45 0.4
Y co-ordinate
0.35 0.3 0.25 0.2 0.15 do diag od diag od1 diag od2 diag
0.1 0.05 0 0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
X velocity
Fig. 6. The control velocities obtained using diag element length definition for = 0.00005 in discretization 3 for test case 1 obtained with optimize-then-discretize and discretize-then-optimize approaches and their variations.
Table 1 Values of he and e for the element with maximum ∇ · a at element center for test case 1 Parameter
adv
diag
he
0.01256 0.1434
0.02795 0.3234
e
From Fig. 6, we see that the principal difference between the controls resulting from the discretizethen-optimze and the optimize-then-discretize approaches with diag element length definition arises from the terms containing ∇ · a in the latter. Such terms manifest themselves due to stabilization and therefore the difference in the results varies with the choice of element length. The maximum and minimum values of ∇ · a, sampled at the center of each element in , are 0.8368 and −0.6348. From Fig. 4, one can also note that for the specific choice of the discretize-then-optimize approach with the diag element length definition, the controls do not exhibit significant oscillations. From Table 1, one can see that with the choice of diag element length, the element length he and the stabilization parameter e are both larger in comparison with the adv element length choice. 4.2. Test case 2 This test case involves flow past a circular cylinder [12, p. 19]. A schematic of the geometry is given in Fig. 7. We use the length L = 10.0 and the radius r = L/20. The advective velocity field a is obtained by solving (36) with Reynolds number Re = |U |D/ = 100, where D is the diameter of the cylinder and U is the inflow velocity. The boundary conditions for a are identical to those indicated in Fig. 7 with a = (U, V )T and Tx and Ty being the x and y component of (a, p) · n, respectively.
246
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Fig. 7. Computational domain and boundary conditions for test case 2, with L = 10 and r = L/20.
Fig. 8. Finite element mesh for test case 2, discretization 1.
For the optimal control problem, we use the objective function (5), where ∗ is defined as the entire computational domain, excluding the first layer of elements at the inflow. Suction and blowing control is applied in normal direction at the back of the half-cylinder (top right quadrant of the cylinder). The top node in the discretization of the control is forced to have zero control velocity. The inflow velocity is u(y) = (0.1, 0)T . All other boundary conditions for u are indicated in Fig. 7. We present results for two dynamic viscosity coefficients = 0.0005 and = 0.00005, and for three discretizations. The coarse discretization consists of 2720 triangular elements, and is shown in Fig. 8—we refer to it as discretization 1. There are 20 control degrees of freedom for this discretization. Discretization 2 consists of 10,880 triangular elements, is a refinement of discretization 1, and has 40 control degrees of freedom. Discretization 3 consists of 24,480 triangular elements, and has 60 control degrees of freedom. Similar to the first test case, we observe from Fig. 9 that there is a strong oscillatory behavior in the controls for the od diag case, especially for the smallest viscosity, i.e., = 0.00005. Fig. 10 shows the dependence of control velocities on the factor , with increasing oscillatory behavior as increases. From Fig. 11, we observe that the sensitivity of the computed control with respect to is large when the diag element length definition is chosen. In particular, for the od approach with diag element length definition large sensitivities are observed in the region that correlates with the spikes in the computed control shown in Fig. 9. Again, similar to test case 1, we see from Fig. 12 that the principal difference between the controls resulting from the do and the od approaches with diag element length definition arises from the
0.6
0.25
circumferential distance from bottom of control surface
do od do od
0.7
0.25
circumferential distance from bottom of control surface
Controlling velocity for µ = 0.0005
0.8
0.25
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 0
0.05
0.1 0.15 X velocity
0.2
Controlling velocity for µ = 0.0005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0
0
0.05
0.1 0.15 X velocity
0.2
Controlling velocity for µ = 0.0005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0
0
0.05
0.1 0.15 X velocity
0.2
247
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Fig. 9. The control velocities obtained for the two optimize-then-discretize and discretize-then-optimize approaches using the diag and adv element length strategies for test case 2. The ordering of these plots is = 0.0005 on the left and = 0.00005 on the right with discretizations 1–3 from top to bottom.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Controlling velocity for µ = 0.00005
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.4
od diag (α = 0.5) od diag (α = 1.0) od diag (α = 2.0) od adv (α = 1.0)
−0.3
−0.2
−0.1 0 X velocity
0.1
0.2
0.3
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
248
Controlling velocity for µ = 0.00005 0.8 0.7 0.6 0.5 0.4 0.3
do diag (α = 0.5) do diag (α = 1.0) do diag (α = 2.0) do adv (α = 1.0)
0.2 0.1 0
0
0.05
0.1
0.15 0.2 X velocity
0.25
0.3
0.35
Fig. 10. The control velocities obtained using different choices of using diag element length definition for = 0.00005 in discretization 3 for test case 2. The plot on the left uses optimize-then-discretize approach and the one on right uses discretize-then-optimize approach. For comparison we also present adv element length definition with = 1.
Table 2 Values of he and e for the element with maximum ∇ · a at element center for test case 2 Parameter
adv
diag
he
0.00694 0.0349
0.0496 0.2522
e
terms containing ∇ · a in the latter. The maximum and minimum values of ∇ · a, sampled at the center of each element in , are 1.2745 and −1.126. As noted in the first test case, one can see from Table 2 that with the choice of diag element length, the element length he and the stabilization parameter e are both larger in comparison with the adv element length choice.
5. Conclusions The order in which a stabilized finite element discretization is applied to an optimal control problem can have a significant effect on the computed solution. If stabilized finite elements are used, the discretizethen-optimize and the optimize-then-discretize approach lead to different optimality conditions on a fixed grid. These differences primarily occur in the discrete adjoint equations and are due to the way the stabilization enters the formulation. Depending on the choice of the element length used for evaluating the stabilization parameter, the computed control may be very sensitive to the choice of the weighting parameter used in the stabilization term and it may exhibit large spurious oscillations. It is important to understand and diagnose these effects. Otherwise, it is questionable whether the computed control resembles the true control.
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5 2 2.5 Weighted sensitivity
3
3.5
circumferential distance from bottom of control surface
Sensitivity plot µ = 0.00005, α = 0.5 0.8
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
249
Sensitivity plot µ = 0.00005, α = 1.0 0.8
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
100
200
300 400 500 Weighted sensitivity
600
700
Sensitivity plot µ = 0.00005, α = 2.0 0.8
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
50
100 150 200 Weighted sensitivity
250
300
Fig. 11. The weighted sensitivity |(g ( ))i |/|g( )i | of the controls for the discretize-then-optimize and optimize-then-discretize approaches using the diag and adv element length definitions for = 0.00005 and discretization 3 for test case 2, with = 0.5, 1.0 and 2.0.
Our computations indicate that a choice (22) of the element length based on the direction of the advective field leads to computed controls that appear to converge as the grid is refined and that, for sufficiently fine grids, are almost independent of the order in which the stabilized finite element discretization is applied to the optimal control. In addition, the sensitivity of the computed control with respect to the weighting parameter used in the stabilization term is small. Our numerical results suggest two indicators for the quality of the computed solution. First, the difference between the controls computed by the discretize-then-optimize approach and the optimize-thendiscretize approach should be small. Secondly, the sensitivity of the computed control with respect to the stabilization parameter should be small. If the discretize-then-optimize approach is chosen as the
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 circumferential distance from bottom of control surface
250
Controlling velocity for µ = 0.00005 0.8
do diag od diag od1 diag od2 diag
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Fig. 12. The control velocities obtained using diag element length definition for = 0.00005 in discretization 3 for test case 2, on experimenting with the terms of difference between the optimal-then-discretize and discretize-then-optimize approaches.
principal solution approach, the system corresponding to the optimize-then-discretize approach and the sensitivity equations can be obtained by minor modifications. Acknowledgements The authors gratefully acknowledge computing resources made available by the National Partnership for Advanced Computational Infrastructure (NPACI). Additional computing resources were provided by the NSF MRI award EIA-0116289. This work was supported by the National Science Foundation under award ACI-0121360, CTS-ITR-0312764, and by Texas ATP Grant 003604-0011-2001. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 32 (1982) 199–259. [2] T.J.R. Hughes, L.P. Franca, M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Meth. Appl. Mech. Eng 59 (1986) 85–99. [3] F. Shakib, Finite element analysis of the compressible Euler and Navier–Stokes equations, Ph.D. Thesis, Department of Mechanical Engineering, Stanford University, 1988. [4] L.P. Franca, S.L. Frey, Stabilized finite element methods: II. The incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 99 (1992) 209–233. [5] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Meth. Appl. Mech. Eng. 73 (1989) 173–189. [6] M. Heinkenschloss, Formulation and analysis of a sequential quadratic programming method for the optimal boundary control of Navier–Stokes flow, in: W.W. Hager, P.M. Pardalos (Eds.), Optimal Control: Theory, Algorithms, and Applications, Kluwer Academic Publishers BV., Dordrecht, 1998, pp. 178–203.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
251
[7] L.S. Hou, S.S. Ravindran, Numerical approximation of optimal control problems by a penalty method: error estimates and numerical results, SIAM J. Sci. Comput. 20 (1999) 1753–1777. [8] S.S. Collis, M. Heinkenschloss, Analysis of the streamline-upwind/Petrov–Galerkin method applied to the solution of optimal control problems’, Technical Report TR02-01, Department of Computational and Applied Mathematics, Rice University, Houston, TX 77005-1892, 2002, http://www.caam.rice.edu/∼heinken/papers/ supg_analysis.html. [9] L.B. Leverett, J.D. Hellums, C.P. Alfrey, E.C. Lynch, Red blood cell damage by shear stress, Biophys. J. 12 (1972) 257–273. [10] M.H. Kroll, J.D. Hellums, L.V. McIntire, A.I. Schafer, J.L. Moake, Platelets and shear stress, Blood 85 (1996) 1525–1541. [11] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [12] O. Ghattas, J.H. Bark, Optimal control of two- and three-dimensional Navier–Stokes flow, J. Computat. Phys. 136 (1997) 231–244. [13] M.D. Gunzburger, L.S. Hou, Treating inhomogeneous essential boundary conditions in finite element methods and the calculation of boundary stresses, SIAM J. Numerical Anal. 29 (1992) 390–424. [14] L.S. Hou, S.S. Ravindran, A penalized Neumann control approach for solving an optimal Dirichlet control problem for the for the Navier–Stokes equations, SIAM J. Control Optim. 36 (1998) 1795–1814. [15] V. Girault, P.A. Raviart, Finite Element Methods for the Navier–Stokes Equations, Springer, New York, 1986. [16] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Computational Mathematics, vol. 15, Springer, New York, 1991. [17] T.E. Tezduyar, Y.J. Park, Discontinuity capturing finite element formulations for nonlinear convection-diffusion-reaction problems, Comput. Meth. Appl. Mech. Eng. 59 (1986) 307–325. [18] K.E. Jansen, S.S. Collis, C. Whiting, F. Shakib, A better consistency for low-order stabilized finite element methods, Comput. Meth. Appl. Mech. Eng. 174 (1999) 153–170.
The effect of stabilization in finite element methods for the optimal boundary control of the Oseen equations Feby Abrahama , Marek Behrb,∗ , Matthias Heinkenschlossc a Department of Mechanical Engineering and Materials Science, Rice University MS 321, 6100 Main Street, Houston,
TX 77005, USA b Lehrstuhl für Numerische Mechanik, Technische Universität München, Boltzmannstr. 15, D-85747 Garching, Germany c Department of Computational and Applied Mathematics, Rice University MS 134, 6100 Main Street, Houston, TX 77005, USA
Received 26 May 2003; received in revised form 14 May 2004; accepted 15 June 2004
Abstract We study the effect of the Galerkin/Least-Squares (GLS) stabilization on the finite element discretization of optimal control problems governed by the linear Oseen equations. Control is applied in the form of suction or blowing on part of the boundary. Two ways of including the GLS stabilization into the discretization of the optimal control problem are discussed. In one case the optimal control problem is first discretized and the resulting finitedimensional problem is then solved. In the other case, the optimality conditions are first formulated on the differential equation level and are then discretized. Both approaches lead to different discrete adjoint equations and, depending on the choice of the stabilization parameters and grid size, may significantly affect the computed control. The effect of the order in which the discretization is applied and the choice of the stabilization parameters are illustrated using two test problems. The cause of the differences in the computed controls are explored numerically. Diagnostics are introduced that may guide the selection of sensible stabilization parameters. 䉷 2004 Elsevier B.V. All rights reserved. Keywords: Optimal boundary control; Stabilized finite element methods; Oseen equations; Solution accuracy
1. Introduction Stabilized finite element methods (FEMs) are frequently and successfully used to discretize advection-dominated partial differential equations (PDEs) [1], or to circumvent the compatibility ∗ Corresponding author.
E-mail addresses: [email protected] (F. Abraham), [email protected] (M. Behr), [email protected] (M. Heinkenschloss). 0168-874X/$ - see front matter 䉷 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2004.06.001
230
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
conditions restricting the choice of interpolation function spaces [2]. One question that arises in the application of these methods is the choice of the stabilization parameter. For many PDE model problems this issue has been studied analytically and numerically [3,4]. However, the impact on the choice of the stabilization parameter when the stabilized FEM is used in the context of optimal control is not well studied. It is not true that a scheme which gives good approximations to the PDE solution for a fixed simulation is also guaranteed to provide good approximations in the context of optimal control problems. Solutions of optimal control problems are characterized by the original governing PDE as well as another PDE—the so-called adjoint equation. Depending on the approach chosen, a discretization of the original governing PDE may imply a discretization scheme for the adjoint equation with poor approximation properties. The present paper investigates this issue for a class of boundary control problems discretized using the Galerkin/least-squares (GLS) approach [5]. The purpose of this paper is to numerically investigate the influence of stabilization parameters on the solution of optimal control problems related to the optimal boundary control of the Navier–Stokes equations. We will demonstrate that the way the stabilization is introduced into the formulation of the optimal control problem, and the choice of stabilization parameters, can have a significant impact on the computed control. Furthermore, we point out some diagnostic tools that may help to assess the quality of the computed control and guide the choice of stabilization parameters. We consider a class of linear quadratic optimal control problems governed by the Oseen equations, and we study the effect of the GLS-stabilized finite element method on the computed control. The linear Oseen equations were chosen as the governing state equation instead of the nonlinear Navier–Stokes equations, because the resulting optimal control problem has a unique solution and the first-order optimality conditions are necessary and sufficient. Optimal control problems governed by the nonlinear Navier–Stokes equations may have local solutions and their solution requires iterative methods. Since it is difficult to separate the possible effects of local solutions and iterative solvers on the computed optimal control from the effects of the stabilization on the computed optimal control, we have chosen the Oseen equations. However, the linear quadratic optimal control problems studied in this paper are closely related to the subproblems that arise in the solution of boundary control problems governed by the Navier–Stokes equations using Newton or sequential quadratic programming methods (see, e.g., [6,7]). Consequently, the results reported in this paper are also relevant for the optimal control of Navier–Stokes flow. A related paper [8] studies the effect of the streamline-upwind/Petrov–Galerkin stabilized finite element method on the numerical solution of linear quadratic distributed optimal control problems governed by an advection-diffusion equation. That paper contains both analytical results that describe the convergence of the computed control (state/adjoint) to the exact control (state/adjoint) as the grid is refined, as well as numerical convergence studies for the simpler optimal control problem. The analytical results in [8] are accurate asymptotically, but do not describe the numerical behavior well for first-order finite elements on relatively coarse, but for practical purposes acceptable grids. Therefore, the present paper focuses on a numerical study. We expect that most of the analytical results in [8] can be extended to the problems and the GLS stabilization method considered in this paper. However, such a theoretical treatment is beyond the scope of this paper. We consider a viscous incompressible fluid occupying a bounded region ⊂ Rnsd , where nsd is the number of space dimensions. The symbols u and p represent the velocity and pressure. The boundary * of is decomposed into three disjoint segments h , d and g . Suction and blowing control is applied
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
on g . We are interested in the solution of the following problem: minimize 2 J (u, g) = J1 (u) + |g|2 dx 2 g
231
(1)
subject to (a · ∇)u − ∇ · (u, p) = f ∇ ·u=0
in ,
(2a)
in ,
(u, p) · n = h
(2b) on h ,
(2c)
u=d
on d ,
(2d)
u=g
on g .
(2e)
The above equations represent the momentum and continuity equations subject to Neumann- and Dirichlettype boundary conditions. In (2), the functions f, h and d are given, and the control g has to be determined as the solution of the optimization problem. The regularization parameter 2 > 0 is given. The stress tensor is given by (u, p) = −pI + 2(u),
(u) =
1 (∇u + ∇uT ), 2
(3)
where is dynamic viscosity and I denotes the identity tensor. We consider two objective functions J (u, g) with 1 J1 (u) = |∇ × u|2 dx (4) 2 ∗ and
J1 (u) = 2
∗
(u): (u) dx,
(5)
where ∗ ⊂ . These objective functions—vorticity magnitude and dissipation due to viscous stress—represent design goals encountered in biomedical engineering. For example, artificial flow devices such as blood pumps induce non-physiological levels of shear stress, leading to blood damage [9] or thrombus formation (clotting) [10]. We give a precise statement of the optimal control problem in Section 2.1, as well as establish existence of a unique solution and formulate the first-order necessary and sufficient optimality conditions. In Section 2.2 we recall the GLS-stabilized finite element method for the discretization of the state equation. The need for stabilization in (2) is twofold. Firstly, GLS removes the under-diffusivity of Galerkin methods which leads to oscillatory solutions in advection-dominated flows. Secondly, GLS circumvents the Babuška–Brezzi (inf–sup) condition, and allows the use of equal-order finite elements for velocities and pressure. For the numerical solution of linear quadratic optimal control problems, there are at least two approaches. In the first approach, called discretize-then-optimize (do), the objective function and governing
232
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
PDEs are discretized. In our case, the Oseen equation is discretized using the GLS-stabilized finite element formulation. This leads to a large-scale finite-dimensional quadratic programming problem, or equivalently, to a large-scale linear system that is solved using suitable numerical linear algebra tools. In the second approach, called optimize-then-discretize (od), the first-order necessary and sufficient optimality conditions are formed and then discretized. The first-order optimality conditions consist of the Oseen equations, the so-called adjoint PDEs—which have a structure similar to the Oseen equations—and an algebraic equation that links controls and adjoint variables. These PDEs are then individually discretized, in our case using the GLS formulation applied to the Oseen and the adjoint PDEs. This also results in a large-scale linear system. These two approaches will be described in Sections 2.3 and 2.4, respectively. The linear systems arising in either approach are not the same, due to the way the stabilization terms affect the adjoint PDEs. These differences and the choice of the GLS stabilization parameter will be explored numerically in Section 4. We will see that an inappropriate choice of the stabilization parameter leads to significant differences in the computed controls. Moreover, the computed control can be very sensitive to a scalar weighting parameter in the stabilization term. In our results, the solutions computed by the optimize-then-discretize approach are more sensitive to the choice of the stabilization parameter and the scalar weighting term. If the stabilization is computed using an element length based on the direction of the advective field a, the optimal controls computed by either approach are in good agreement on fine grids. 2. A model problem 2.1. Formulation of the optimal control problem In the following, we use the spaces L2 () and H 1 (), which are defined in the usual way [11], and V = {v ∈ H 1 () | v = 0 on d }. We set H1 () = [H 1 ()]nsd , V = V nsd , L2 (h ) = [L2 (h )]nsd , and L2 (g ) = [L2 (g )]nsd . The Dirichlet boundary conditions (2d) and (2e) can be implemented through interpolation [12], weakly through a Lagrange multiplier technique [13], or via a penalty approach [7,14]. We use interpolation to implement the Dirichlet boundary conditions (2d) with fixed data, and replace (2e) by the penalized Neumann boundary condition (u, p) · n +
1
u=
1
g
on g ,
(6)
where > 0 is a given penalty parameter. This choice enables us to look for controls g in L2 (g ), instead of H1/2 (g ). The weak form of the partial differential equation (2a)–(2d), (6) is given as follows: Find u ∈ H1 () with u = d on d and p ∈ L2 (), such that 1 (a · ∇)u · v dx + (v): (u, p) dx − q∇ · u dx + v · (u − g) dx g = v · h dx + f · v dx (7) h
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
233
for all v ∈ V and all q ∈ L2 (). It is also possible to implement (2d) via a penalized Neumann approach. See Remark 2.3 below. Eq. (7) motivates the definition of the bilinear forms a(u, v) = 2 (v): (u) dx ∀v, u ∈ H1 (), (8a)
c(u, v) =
(a · ∇)u · v dx
∀v, u ∈ H1 (),
(8b)
b(u, q) =
q∇ · u dx
∀u ∈ H1 (), ∀q ∈ L2 (),
(8c)
h · v dx
∀h ∈ L2 (h ), ∀v ∈ V,
(8d)
g · v dx
∀g ∈ L2 (h ), ∀v ∈ V
(8e)
h, vh = g, vg =
h
g
and of the linear functional f, v = f · v dx ∀v ∈ V.
(8f)
For a given > 0, we consider the optimal control problem 2 minimize J1 (u) + |g|2 dx, 2 g subject to
(9a)
a(u, v) + c(u, v) − b(v, p) − b(u, q) + −1 u, vg − −1 g, vg = h, vh + f, v
∀v ∈ V, ∀q ∈ L2 ().
(9b)
Lemma 2.1. Let > 0 be given and assume that there exists ud ∈ H1 () with ud =d on d . If a ∈ H1 (h ) satisfies ∇ · a = 0 and a · n 0 on \d , then for any h ∈ L2 (h ), g ∈ L2 (g ) and f ∈ L2 () Eq. (9b) has a unique solution (u, p) ∈ H1 () × L2 (). Moreover, there exists a constant C > 0 (dependent on , but independent of h, g, f and ud ) such that uH1 () + pL2 () C(hL2 (h ) + gL2 (h ) + fL2 () + ud H1 () ). Proof. The bilinear form b satisfies the inf–sup condition inf
sup b(v, q)/(qL2 () vV ) > 0
q∈L2 () v∈V
[15, p. 81]. Integration by parts yields 1 1 2 c(v, v) = − ∇ · a |v| dx + a · n |v|2 dx 2 2 * With the assumptions on a, this implies c(v, v) 0
∀v ∈ V.
∀v ∈ H1 ().
(10)
234
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Hence, there exists > 0 such that a(v, v) + c(v, v) + −1 v, vg v2H1 ()
∀v ∈ V.
(11)
Continuity of the linear and bilinear forms (8) follows using the arguments in [14]. The assertion now follows from the theory of saddle point problems [16, Section II.1], [15, Section 4]. Theorem 2.2. Let the assumptions of Lemma 2.1 be satisfied. If J1 is weakly lower semicontinuous on V, then (9) admits a solution (u , p , g ) ∈ H1 () × L2 () × L2 (g ). If J1 is convex, e.g., if J1 is given by (4) or (5), the solution is unique. Proof. The result follows from standard arguments, see, e.g., the proof of Theorem 3.5 in [14]. Remark 2.3. It is possible to also implement the Dirichlet boundary condition (2d) via a penalized Neumann condition 1 1 1 (u, p) · n − (a · n) u + u = d on d . (12) 2 The weak form of the partial differential equation (2a)–(2c), (6), (12) is given as follows: Find u ∈ H1 () and p ∈ L2 (), such that 1 (a · ∇)u · v dx + (v): (u, p) dx − q∇ · u dx − (a · n)|u|2 dx 2 d 1 1 v · (u − g) dx + v · (u − d) dx + g d = v · h dx + f · v dx (13) h
for all v ∈ and all q ∈ L2 (). Lemma 2.1 with V replaced by H1 () still holds, since (10) and the inclusion of the term − 21 (a · n) u in (12) imply H1 ()
a(v, v) + c(v, v) + −1 v, vg v2H1 ()
∀v ∈ H1 ().
(14)
The convergence behavior of the solutions (u , p , g ) of the optimal control problem (9) with penalized Neumann boundary control to the solution of (1) and (2) as → 0 is discussed in [14] in the context of Navier–Stokes equations. For the purpose of this study, we view (9) with a small but fixed > 0 as our optimal control problem. The Lagrangian associated with the optimal control problem (9) is given by 2 L(u, p, g, , )=J1 (u) + |g|2 dx 2 g + a(u, ) + c(u, ) − b(, p) − b(u, ) + −1 u, g − −1 g, g − h, h − f, . (15) Lemma 2.1 ensures that (9) satisfies a constraint qualification. Hence, the necessary and sufficient optimality conditions for (9) are obtained by setting the Fréchet derivatives of the Lagrangian with
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
235
respect to the adjoint variables (, ), with respect to the state variables (u, p) and with respect to the controls g equal to zero. The necessary and sufficient optimality conditions consist of the state equation a(u, v) + c(u, v) − b(v, p) − b(u, q) + −1 u, vg − −1 g, vg = h, vh + f, v
∀v ∈ V, ∀q ∈ L2 (),
(16a)
the adjoint equation a(v, ) + c(v, ) − b(, q) − b(v, ) + −1 v, g = −Du J1 (u), v
∀v ∈ V, ∀q ∈ L2 ()
(16b)
∀z ∈ L2 (g ).
(16c)
and the gradient equation 2 g, zg =
1
, zg
In the adjoint equation (16b), Du J1 (u) denotes the Fréchet derivative of J1 with respect to u and ·, · is the duality pairing between (H1 ())∗ —the dual of H1 ()—and H1 (), i.e., Du J1 (u), v is the application of the Fréchet derivative of J1 to the function v. In the case of our objective functions (4) and (5) these are simply given as Du J1 (u), v = (∇ × u) · (∇ × v) dx (17) ∗
and
Du J1 (u), v = 4
∗
(u): (v) dx,
(18)
respectively. The adjoint equation (16b) may formally be interpreted as the weak form of the adjoint partial differential equation −(a · ∇) − (∇ · a) − ∇ · (, ) = −Du J1 (u) ∇ ·=0
in ,
1 2
on g .
(19c) (19d)
1
and Eq. (16c) simply states that g=
on h ,
on d ,
(, ) · n + (a · n) +
(19a) (19b)
(, ) · n + (a · n) = 0 =0
in ,
=0
on g ,
(19e)
236
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
2.2. Discretization of the state equations For the discretization of the state equation we apply the GLS-stabilized finite element method using piece-wise-linear polynomials for both the velocity and the pressure. We divide our domain into nel subdomains e , e = 1, 2, . . . , nel . We assume that our triangulation is such that the controlled boundary g can be written as the union of edges or faces of elements e , e = 1, 2, . . . , nel . We use piecewiselinear polynomials to approximate both the velocities and pressures. Hence, the discretized velocities are in H1h () = [H 1h ()]nsd , where H 1h () = { h ∈ C 0 () | h |e ∈ P 1 , e = 1, 2, . . . , nel } and the discretized pressures are in L2h () = { h ∈ C 0 () | h |e ∈ P 1 , e = 1, 2, . . . , nel }. Furthermore, we define V h = { h ∈ H 1h () | h = 0 on d }, Vh = [V h ]nsd , L2h (g ) = { h |g | h ∈ H 1h ()} and L2h (g ) = [L2h (g )]nsd . We let Ih d denote the piecewise linear interpolant of the Dirichlet boundary data d. The GLS-stabilized finite element formulation for the state equation (2a)–(2c), (6) is given as follows: Find uh ∈ H1h () with uh = Ih d on d and ph ∈ L2h () such that h h v · (a · ∇)u dx + (vh ) : (uh , ph ) dx
+
nel e
e [(a · ∇)vh − ∇ · (vh , q h )] · [(a · ∇)uh − ∇ · (uh , p h )] dx
e=1 1 h h + q ∇ · u dx + =
h
vh · h dx +
g
vh · (uh − gh ) dx
f · vh dx +
nel e=1
e
e [(a · ∇)vh − ∇ · (vh , q h )] · f dx
(20)
for all vh ∈ Vh and all q h ∈ L2h (). Our choice [3] for the stabilization parameter e is
e =
2|a| he
2
+
4 h2e
2 −1/2 ,
(21)
where is a factor typically taken to be 1.0 and the element length he is computed using one of two choices diag and adv . Let h∗e be the longest edge in an element. The diag element length is given by he = h∗e . The adv element length, introduced in [17], is defined as ∗ if |aeh | = 0, he
−1 he = (22) otherwise, a |se · ∇Na | where aeh is the interpolated advection velocity at the center of element e, se = aeh /|aeh |, and Na represents the basis functions. Both these element length definitions are commonly used in the computation of
e . Notice that the adv element length is identical to the diag element length in regions where there is stagnation in the advection field. In many flow scenarios with extended and frequent regions of stagnation, the diag element length is used to avoid abrupt switching of the element length type from streamline interpolation to longest edge of the element, as is the case in the adv element length definition.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
237
2.3. Discretization of the optimization problem A frequently used approach for the numerical solution of an optimal control problem is to discretize the optimal control problem first and then solve the resulting finite-dimensional nonlinear (in our case quadratic) programming problem using a suitable optimization algorithm. We refer to this approach as the discretize-then-optimize (do) approach. Using the GLS-stabilized finite element method described in Section 2.2 for the discretization of the state equation, our discretization of the optimal control problem (9) is given by 2 J1 (uh ) + |gh |2 dx, 2 g (20), uh ∈ H1h (), ph ∈ L2h (), gh ∈ L2h (g ), uh = Ih d on d .
minimize subject to with
(23)
The Lagrangian for the discretized problem (23) is given by Lh (uh , ph , gh , h , h ) 2 h = J1 (u ) + |gh |2 dx 2 g h h + · (a · ∇)u dx + (h ): (uh , ph ) dx
+
e=1 + −
nel
e
h
e [(a · ∇)h − ∇ · (h , h )] · [(a · ∇)uh − ∇ · (uh , p h )] dx h
∇ · u dx +
h
h · h dx −
1
g
h · (uh − gh ) dx
f · h dx −
nel e=1
e
e [(a · ∇)h − ∇ · (h , h )] · f dx.
(24)
The necessary and sufficient optimality conditions for (23) are obtained by setting the derivatives of the Lagrangian with respect to the discrete adjoint variables (h , h ), the discretized state variables (uh , ph ) and the discretized controls gh to zero. Setting the derivative of Lh with respect to the discrete adjoint variables (h , h ) to zero gives the discretized state equation (20). Setting the derivative of Lh with respect to the discretized state variables (uh , ph ) to zero gives the discrete adjoint equation h h · (a · ∇)v dx + (h ): (vh , q h ) dx
+
e=1 +
nel
e
h
e [(a · ∇)h − ∇ · (h , h )] · [(a · ∇)vh − ∇ · (vh , q h )] dx h
∇ · v dx +
= −Du J1 (uh ), vh
1
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(25)
238
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
where Du J1 (uh ), vh is defined as in (17) or (18), respectively. Finally, setting the derivative of Lh with respect to the discretized control variables gh to zero gives the discrete gradient equation 1 2 h h g · z dx = h · zh dx ∀zh ∈ L2h (g ). (26)
g
g
The necessary and sufficient conditions (25), (26) and (20) lead to a system of linear equations
Hdo 0 A
0 GT −B
AT 0 u do T gdo = 0 . −B f 0 do
(27)
Here u do and do are the vectors containing the values of (uh , ph ) and (h , h ) at the grid points, respectively. 2.4. Discretization of the optimality conditions The optimality conditions (16) are necessary and sufficient. Therefore, we can tackle the optimality conditions directly to compute an approximate solution of our optimal control problem. This is called the optimize-then-discretize (od) approach. This approach requires us to discretize the state and the adjoint PDE (19). In principle, these discretizations do not have to be the same. However, we will use the same mesh and the GLS-stabilized finite element method for the discretization of both. We use the same notation as in Section 2.2. Our discretization of the state equation is given in (20). The GLS-stabilized finite element method applied to the adjoint equation (19), yields h h · (a · ∇)v dx + (vh ): (h , h ) dx
+
nel e
e [−(a · ∇)vh − (∇ · a)vh − ∇ · (vh , q h )] · [−(a · ∇)h − (∇ · a)h − ∇ · (h , h )] dx
e=1 1 h h + q ∇ · dx +
= −Du J1 (uh ), vh s
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(28)
where Du J1 (uh ), vh s is given as follows. For the objective function (4), Du J1 (uh ), vh s = (∇ × uh ) · (∇ × vh ) dx ∗ nel
+
e=1
e ∩∗
e (∇ × uh ) · (∇ × [−(a · ∇)vh
− (∇ · a)vh − ∇ · (vh , q h )]) dx,
(29)
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
239
and for the objective function (5), Du J1 (uh ), vh s =4 (uh ): (vh ) dx ∗
+ 4
nel e=1
e ∩∗
e (uh ): (−(a · ∇)vh
− (∇ · a)vh − ∇ · (vh , q h )) dx.
(30)
The discretization of (16c) is again given by (26). The necessary and sufficient conditions (20), (28) and (26) lead to a system of linear equations 0 A 0 Hod u od T T 0 (31) G −B god = 0 . f A −B 0 od Here u od and od are the vectors containing the values of (uh , ph ) and (h , h ) at the grid points, respectively. Note that because of the stabilization terms in (29) and in (30), the submatrix Hdo in (27) is different from the submatrix Hod in (31). Moreover, because the advective terms in adjoint equation (19) are different from those in the Oseen equation (2), the discrete adjoint equations (25) are different from the discretized adjoint equations (28). Hence the submatrix AT in (27) is also different from submatrix A in (31). 3. Implementation 3.1. Stabilization terms Our choice of piecewise-linear functions makes our discrete approximation a low-order method, i.e. the order of the function space is lower than the order of the highest derivative in (2), the stabilized finite element equivalent of which is (20). Therefore, the terms ∇ · (vh , q h ) and ∇ · (uh , ph ) in (20) are given by ∇q h and ∇p h , respectively. Dropping the term ∇ · 2(uh ) guarantees us only a weak consistency [18], since this term itself vanishes with refinement, as e approaches zero based on the relation in (21). Normally, a variational reconstruction of this term is done to achieve a “stronger consistency” for such lower-order choice of function space [18]. However, we have not done so here, to avoid the potential nonlinearity thus introduced. A similar inconsistency is introduced into (25) and (28). The implemented version of (25) is given by h h · (a · ∇)v dx + (h ): (vh , q h ) dx
+
e=1 +
nel
e
h
e [(a · ∇)vh + ∇q h ] · [(a · ∇)h + ∇ h ] dx h
∇ · v dx +
= −Du J1 (uh ), vh
1
g h
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
(32)
240
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
and the implemented version of (28) is given by h h · (a · ∇)v dx + (vh ): (h , h ) dx
+
nel e
e [−(a · ∇)vh − (∇ · a)vh + ∇q h ] · [−(a · ∇)h − (∇ · a)h + ∇ h ] dx
e=1 1 h h + q ∇ · dx +
= −Du J1 (uh ), vh s where h
h
Du J1 (u ), v s =
∗
+
h
h
(33)
(∇ × uh ) · (∇ × vh ) dx
nel
Du J1 (u ), v s =4
h · vh dx
∀v ∈ Vh , ∀q ∈ L2h (),
e=1
and
g h
∗
+ 4
e ∩∗
e (∇ × uh ) · (∇ × [−(a · ∇)vh − (∇ · a)vh ]) dx,
(34)
(uh ): (vh ) dx
nel e=1
e ∩∗
e (uh ): (−(a · ∇)vh − (∇ · a)vh ) dx,
(35)
for the objective functions (4) and (5), respectively. 3.2. Differences The use of the GLS-stabilized finite element method in the discretization of the optimal control problem creates differences between the discretize-then-optimize approach and the optimize-then-discretize approach. Specifically, there are differences between the discrete adjoint equation (32) and the discretized adjoint equation (33). To explore how these differences impact the computed solution, we also implement variations of (33). 1. The left- and right-hand side in (33) contain ∇ · a terms, while the divergence of the advective field does not enter into the discrete adjoint equation (32). Therefore, we also compute the optimal controls using (20), (26) and (33), where the ∇ · a terms in (33) are dropped. The optimal controls computed this way will be labeled as od1. 2. If we compare (32) and (33) with ∇ · a terms in (33) replaced by zero, we see that the right-hand sides of the discretized adjoint equations (33) contain −(a · ∇)vh terms that are not present in the discrete adjoint equations (32). Therefore, we also compute the optimal controls using (20), (26) and (33) where all ∇ · a terms in (33) as well as the −(a · ∇)vh term on the right-hand side of (33) are dropped. The optimal controls computed this way will be labeled as od2.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
241
3. After the modifications of the discretized adjoint (33) described above in 2. have been performed, the only remaining difference between (32) and (33) is the sign of (a · ∇)vh and (a · ∇)h .
3.3. Sensitivity analysis We are interested in the sensitivity of the computed solution with respect to the stabilization factor . Eqs. (25), (26) and (20) reveal that the submatrix A and the right hand side component f in (27) depend on . Hence we write A( ) and f( ) instead of A and f, respectively. Moreover, we denote the solution of (27) by u do ( ), gdo ( ), do ( ). Examination of (25), (26) and (20) reveals that the submatrices A( ) and Hdo ( ) in (27) are of the form A( ) = A1 + A2 and Hdo ( ) = H1 , where we use subscript 1 for all those terms that do not depend on stabilization, and subscript 2 for those terms that are dependent on ( )=0. stabilization. Hence, if we use to denote differentiation with respect to , then A ( )=A2 and Hdo By the implicit function theorem, the derivatives of the solution of(27) satisfy
Hdo
0
A( )T
0 GT
A( ) −B
( ) u do
−BT gdo ( ) 0 do ( )
0
= − 0 A ( )
0
(A ( ))T
0
0
0
0
0 gdo ( ) + 0 , f ( ) do ( ) u do ( )
where u do ( ), gdo ( ) and do ( ) solve (27). The sensitivity of the solution u od ( ), god ( ) and od ( ) of (31) can be computed analogously. Note that in this case Hod also depends on and is of the form Hod ( ) = H1 + H2 , with Hod ( ) = H2 .
4. Numerical results In this section we study the effect of the GLS stabilization on the computed control for two test cases derived from commonly used model problems. In both cases, the advective field a is computed by solving the Navier–Stokes equations (a · ∇)a − ∇ · (a, p) = 0
∇ · a=0
on ,
on , (36)
with appropriate boundary conditions. The Navier–Stokes equations (36) are discretized using the GLSstabilized finite element method and the resulting nonlinear system is solved using Newton’s method with a residual tolerance of 10−16 . Note that while the “exact” advective field a is divergence-free, this is not true for the computed advective field, which is used as the coefficient in our computations. We solve (9) with parameters = 10−5 and 2 = 10−5 . As stated before, we use the notation do to refer to the control computed using the discretize-then-optimize approach (26), (20) and (32), and od, to refer to the control computed by the optimize-then-discretize approach (26), (20) and (33), respectively.
242
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Γt Γin Γg
Ω
Ω*
Γo
Γb
Fig. 1. Computational domain for test case 1.
Fig. 2. Finite element mesh for test case 1, discretization 1.
4.1. Test case 1 The first problem is modeled after the backward-facing step problem [7, pp. 1767,1769]. A schematic of the geometry is given in Fig. 1. The height of the inflow boundary is 0.5 and that of the outflow boundary is 1. The length of the narrower section of the channel is 1 and that of the wider section is 7, with the total horizontal length being 8. The advective velocity field a is computed by solving the Navier–Stokes equations (36) with the following boundary conditions. The inflow velocity is assumed parabolic with a profile a(0, y) = (8(y − 0.5)(1 − y), 0)T . At the outflow boundary o we impose traction-free boundary conditions for a in the x-direction. No-slip conditions are imposed at all other boundaries. We define the Reynolds number Re = Umax H /, where Umax is the maximum inlet velocity and H is the channel height. The Reynolds number for computing a is Re = 200. The boundary conditions for the Oseen equation are as follows. The inflow velocity field, different from that of a, is taken as parabolic given by u(0, y) = ((y − 0.5)(1 − y), 0)T . A traction-free boundary condition for the x-velocity is imposed at o . The controlled boundary g is the backward facing step. We allow suction and blowing in the normal direction. At the remaining boundaries b and t no-slip conditions are imposed for u. The objective function is given by (4), where ∗ is the subdomain of length 2 and height 0.5, as indicated in Fig. 1. We present results for two dynamic viscosity coefficients = 0.0005 and 0.00005, and for three discretizations. Discretization 1 consists of 1800 triangular elements as shown in Fig. 2. Discretization 2 consisting of 7200 triangular elements and is a uniform refinement of discretization 1. Discretization 3 is an uniform refinement of discretization 2 and consists of 28,800 elements. Fig. 3 shows the control velocities obtained using the various approaches considered here. We observe that the od and do approaches with the adv element length definition converge with mesh refinement. When we use the diag element length, we observe that the discretize-then-optimize approach converges to the solution obtained using the adv element length. However, we can see that there is a strong oscillatory behavior in the controls for optimize-then-discretize diag case, especially for the smallest viscosity, i.e. = 0.00005. Fig. 4 shows the dependence of control velocities on the factor , with increasing oscillatory behavior as increases. Fig. 5 shows the sensitivity of the controls to . The sensitivity of the computed control with respect to is large when the diag element length definition is chosen. In particular for the od approach with diag element length definition large sensitivities are observed near the corner of
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Controlling velocity for µ = 0.0005
0.5
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
Controlling velocity for µ = 0.00005
0.5
0.45
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
do od do od
0.1 0.05
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015 0.02 X velocity
0.025
0.03
0.035
0
0.005 0.01
0.5
0.45
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
0.5
0.3 0.25 0.2 0.15 do od do od
0.05
0.035 0.04 0.045
0.3 0.25 0.2 0.15
0.1
0.015 0.02 0.025 0.03 X velocity
Controlling velocity for µ = 0.00005
Controlling velocity for µ = 0.0005
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015 0.02 X velocity
0.025
0.03
0.035
0
Controlling velocity for µ = 0.0005
0.5
0.005
0.01
0.015
0.02 0.025 X velocity
0.03
0.035
0.04
0.045
0.04
0.045
Controlling velocity for µ = 0.00005
0.5
0.45
0.45
0.4
0.4
0.35
0.35 Y co ordinate
Y co ordinate
243
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
do od do od
0.1 0.05
diag diag adv adv
do od do od
0.1
diag diag adv adv
0.05
0
0 0
0.005
0.01
0.015
0.02
X velocity
0.025
0.03
0.035
0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
X velocity
Fig. 3. The control velocities obtained for the two optimize-then-discretize and discretize-then-optimize approaches using the diag and adv element length definition with = 1 for test case 1. The ordering of these plots is = 0.0005 on the left and = 0.00005 on the right with discretizations 1–3 from top to bottom.
the backward facing step. This correlates with the ‘spikes’ in the computed control for this case shown in Fig. 3. The size of the sensitivities are relatively small for both, the discretize-then-optimize and the optimize-then-discretize approach when the adv element length definition is chosen. Again this correlates with the observed grid convergence for both approaches in case of the adv element length definition (see Fig. 3).
244
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Controlling velocity for µ = 0.00005
0.5
0.4
0.35
0.35 Y co-ordinate
Y co-ordinate
0.45
0.4
0.3 0.25 0.2 0.15 0.1 0.05
Controlling velocity for µ = 0.00005
0.5
0.45
od od od od
0.3 0.25 0.2 0.15
diag (α = 0.5) diag (α = 1.0) diag (α = 2.0) adv (α = 1.0)
do do do do
0.1 0.05
0 -0.1 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08 X velocity
diag (α = 0.5) diag (α = 1.0) diag (α = 2.0) adv (α = 1.0)
0 0.1
0
0.005
0.01
0.015
0.02 0.025 X velocity
0.03
0.035
0.04
Fig. 4. The control velocities obtained using different choices of using diag element length definition for = 0.00005 in discretization 3 for test case 1. The plot on the left uses optimize-then-discretize approach and the one on right uses discretize-then-optimize approach. For comparison we also present adv element length definition with = 1.
Sensitivity plot for µ = 0.00005, α = 0.5
0.5 0.45
0.45
0.4
0.4
0.35
0.35 Y co-ordinate
Y co-ordinate
Sensitivity plot for µ = 0.00005, α = 1.0
0.5
0.3 0.25 0.2 0.15
0.3 0.25 0.2 0.15
0.1
0.1
do diag od diag do adv od adv
0.05
do diag od diag do adv od adv
0.05
0
0 0
0.1
0.2
0.3 0.4 0.5 Weighted sensitivity
0
0.7
0.6
0.2
0.4
0.6 0.8 Weighted sensitivity
1
1.2
1.4
Sensitivity plot for µ = 0.00005, α = 2.0
0.5 0.45 0.4
Y co-ordinate
0.35 0.3 0.25 0.2 0.15 0.1
do diag od diag do adv od adv
0.05 0 0
1
2
3
4
5
6
Weighted sensitivity
Fig. 5. The weighted sensitivity |(g ( ))i |/|g( )i | of the controls for the discretize-then-optimize and optimize-then-discretize approaches using the diag and adv element length definitions for = 0.00005 and discretization 3 for test case 1, with = 0.5, = 1.0 and = 2.0.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
245
Controlling velocity for µ = 0.00005 0.5 0.45 0.4
Y co-ordinate
0.35 0.3 0.25 0.2 0.15 do diag od diag od1 diag od2 diag
0.1 0.05 0 0
0.005
0.01
0.015
0.02
0.025
0.03
0.035
0.04
0.045
X velocity
Fig. 6. The control velocities obtained using diag element length definition for = 0.00005 in discretization 3 for test case 1 obtained with optimize-then-discretize and discretize-then-optimize approaches and their variations.
Table 1 Values of he and e for the element with maximum ∇ · a at element center for test case 1 Parameter
adv
diag
he
0.01256 0.1434
0.02795 0.3234
e
From Fig. 6, we see that the principal difference between the controls resulting from the discretizethen-optimze and the optimize-then-discretize approaches with diag element length definition arises from the terms containing ∇ · a in the latter. Such terms manifest themselves due to stabilization and therefore the difference in the results varies with the choice of element length. The maximum and minimum values of ∇ · a, sampled at the center of each element in , are 0.8368 and −0.6348. From Fig. 4, one can also note that for the specific choice of the discretize-then-optimize approach with the diag element length definition, the controls do not exhibit significant oscillations. From Table 1, one can see that with the choice of diag element length, the element length he and the stabilization parameter e are both larger in comparison with the adv element length choice. 4.2. Test case 2 This test case involves flow past a circular cylinder [12, p. 19]. A schematic of the geometry is given in Fig. 7. We use the length L = 10.0 and the radius r = L/20. The advective velocity field a is obtained by solving (36) with Reynolds number Re = |U |D/ = 100, where D is the diameter of the cylinder and U is the inflow velocity. The boundary conditions for a are identical to those indicated in Fig. 7 with a = (U, V )T and Tx and Ty being the x and y component of (a, p) · n, respectively.
246
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
Fig. 7. Computational domain and boundary conditions for test case 2, with L = 10 and r = L/20.
Fig. 8. Finite element mesh for test case 2, discretization 1.
For the optimal control problem, we use the objective function (5), where ∗ is defined as the entire computational domain, excluding the first layer of elements at the inflow. Suction and blowing control is applied in normal direction at the back of the half-cylinder (top right quadrant of the cylinder). The top node in the discretization of the control is forced to have zero control velocity. The inflow velocity is u(y) = (0.1, 0)T . All other boundary conditions for u are indicated in Fig. 7. We present results for two dynamic viscosity coefficients = 0.0005 and = 0.00005, and for three discretizations. The coarse discretization consists of 2720 triangular elements, and is shown in Fig. 8—we refer to it as discretization 1. There are 20 control degrees of freedom for this discretization. Discretization 2 consists of 10,880 triangular elements, is a refinement of discretization 1, and has 40 control degrees of freedom. Discretization 3 consists of 24,480 triangular elements, and has 60 control degrees of freedom. Similar to the first test case, we observe from Fig. 9 that there is a strong oscillatory behavior in the controls for the od diag case, especially for the smallest viscosity, i.e., = 0.00005. Fig. 10 shows the dependence of control velocities on the factor , with increasing oscillatory behavior as increases. From Fig. 11, we observe that the sensitivity of the computed control with respect to is large when the diag element length definition is chosen. In particular, for the od approach with diag element length definition large sensitivities are observed in the region that correlates with the spikes in the computed control shown in Fig. 9. Again, similar to test case 1, we see from Fig. 12 that the principal difference between the controls resulting from the do and the od approaches with diag element length definition arises from the
0.6
0.25
circumferential distance from bottom of control surface
do od do od
0.7
0.25
circumferential distance from bottom of control surface
Controlling velocity for µ = 0.0005
0.8
0.25
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 0
0.05
0.1 0.15 X velocity
0.2
Controlling velocity for µ = 0.0005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0
0
0.05
0.1 0.15 X velocity
0.2
Controlling velocity for µ = 0.0005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0
0
0.05
0.1 0.15 X velocity
0.2
247
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Controlling velocity for µ = 0.00005
0.8
do od do od
0.7 0.6
diag diag adv adv
0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Fig. 9. The control velocities obtained for the two optimize-then-discretize and discretize-then-optimize approaches using the diag and adv element length strategies for test case 2. The ordering of these plots is = 0.0005 on the left and = 0.00005 on the right with discretizations 1–3 from top to bottom.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 Controlling velocity for µ = 0.00005
0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.4
od diag (α = 0.5) od diag (α = 1.0) od diag (α = 2.0) od adv (α = 1.0)
−0.3
−0.2
−0.1 0 X velocity
0.1
0.2
0.3
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
248
Controlling velocity for µ = 0.00005 0.8 0.7 0.6 0.5 0.4 0.3
do diag (α = 0.5) do diag (α = 1.0) do diag (α = 2.0) do adv (α = 1.0)
0.2 0.1 0
0
0.05
0.1
0.15 0.2 X velocity
0.25
0.3
0.35
Fig. 10. The control velocities obtained using different choices of using diag element length definition for = 0.00005 in discretization 3 for test case 2. The plot on the left uses optimize-then-discretize approach and the one on right uses discretize-then-optimize approach. For comparison we also present adv element length definition with = 1.
Table 2 Values of he and e for the element with maximum ∇ · a at element center for test case 2 Parameter
adv
diag
he
0.00694 0.0349
0.0496 0.2522
e
terms containing ∇ · a in the latter. The maximum and minimum values of ∇ · a, sampled at the center of each element in , are 1.2745 and −1.126. As noted in the first test case, one can see from Table 2 that with the choice of diag element length, the element length he and the stabilization parameter e are both larger in comparison with the adv element length choice.
5. Conclusions The order in which a stabilized finite element discretization is applied to an optimal control problem can have a significant effect on the computed solution. If stabilized finite elements are used, the discretizethen-optimize and the optimize-then-discretize approach lead to different optimality conditions on a fixed grid. These differences primarily occur in the discrete adjoint equations and are due to the way the stabilization enters the formulation. Depending on the choice of the element length used for evaluating the stabilization parameter, the computed control may be very sensitive to the choice of the weighting parameter used in the stabilization term and it may exhibit large spurious oscillations. It is important to understand and diagnose these effects. Otherwise, it is questionable whether the computed control resembles the true control.
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.5
1
1.5 2 2.5 Weighted sensitivity
3
3.5
circumferential distance from bottom of control surface
Sensitivity plot µ = 0.00005, α = 0.5 0.8
circumferential distance from bottom of control surface
circumferential distance from bottom of control surface
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
249
Sensitivity plot µ = 0.00005, α = 1.0 0.8
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
100
200
300 400 500 Weighted sensitivity
600
700
Sensitivity plot µ = 0.00005, α = 2.0 0.8
do diag od diag do adv od adv
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
50
100 150 200 Weighted sensitivity
250
300
Fig. 11. The weighted sensitivity |(g ( ))i |/|g( )i | of the controls for the discretize-then-optimize and optimize-then-discretize approaches using the diag and adv element length definitions for = 0.00005 and discretization 3 for test case 2, with = 0.5, 1.0 and 2.0.
Our computations indicate that a choice (22) of the element length based on the direction of the advective field leads to computed controls that appear to converge as the grid is refined and that, for sufficiently fine grids, are almost independent of the order in which the stabilized finite element discretization is applied to the optimal control. In addition, the sensitivity of the computed control with respect to the weighting parameter used in the stabilization term is small. Our numerical results suggest two indicators for the quality of the computed solution. First, the difference between the controls computed by the discretize-then-optimize approach and the optimize-thendiscretize approach should be small. Secondly, the sensitivity of the computed control with respect to the stabilization parameter should be small. If the discretize-then-optimize approach is chosen as the
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251 circumferential distance from bottom of control surface
250
Controlling velocity for µ = 0.00005 0.8
do diag od diag od1 diag od2 diag
0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 −0.1
−0.05
0
0.05
0.1 0.15 X velocity
0.2
0.25
0.3
Fig. 12. The control velocities obtained using diag element length definition for = 0.00005 in discretization 3 for test case 2, on experimenting with the terms of difference between the optimal-then-discretize and discretize-then-optimize approaches.
principal solution approach, the system corresponding to the optimize-then-discretize approach and the sensitivity equations can be obtained by minor modifications. Acknowledgements The authors gratefully acknowledge computing resources made available by the National Partnership for Advanced Computational Infrastructure (NPACI). Additional computing resources were provided by the NSF MRI award EIA-0116289. This work was supported by the National Science Foundation under award ACI-0121360, CTS-ITR-0312764, and by Texas ATP Grant 003604-0011-2001. References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 32 (1982) 199–259. [2] T.J.R. Hughes, L.P. Franca, M. Balestra, A new finite element formulation for computational fluid dynamics: V. Circumventing the Babuška–Brezzi condition: a stable Petrov–Galerkin formulation of the Stokes problem accommodating equal-order interpolations, Comput. Meth. Appl. Mech. Eng 59 (1986) 85–99. [3] F. Shakib, Finite element analysis of the compressible Euler and Navier–Stokes equations, Ph.D. Thesis, Department of Mechanical Engineering, Stanford University, 1988. [4] L.P. Franca, S.L. Frey, Stabilized finite element methods: II. The incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Eng. 99 (1992) 209–233. [5] T.J.R. Hughes, L.P. Franca, G.M. Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/least-squares method for advective-diffusive equations, Comput. Meth. Appl. Mech. Eng. 73 (1989) 173–189. [6] M. Heinkenschloss, Formulation and analysis of a sequential quadratic programming method for the optimal boundary control of Navier–Stokes flow, in: W.W. Hager, P.M. Pardalos (Eds.), Optimal Control: Theory, Algorithms, and Applications, Kluwer Academic Publishers BV., Dordrecht, 1998, pp. 178–203.
F. Abraham et al. / Finite Elements in Analysis and Design 41 (2004) 229 – 251
251
[7] L.S. Hou, S.S. Ravindran, Numerical approximation of optimal control problems by a penalty method: error estimates and numerical results, SIAM J. Sci. Comput. 20 (1999) 1753–1777. [8] S.S. Collis, M. Heinkenschloss, Analysis of the streamline-upwind/Petrov–Galerkin method applied to the solution of optimal control problems’, Technical Report TR02-01, Department of Computational and Applied Mathematics, Rice University, Houston, TX 77005-1892, 2002, http://www.caam.rice.edu/∼heinken/papers/ supg_analysis.html. [9] L.B. Leverett, J.D. Hellums, C.P. Alfrey, E.C. Lynch, Red blood cell damage by shear stress, Biophys. J. 12 (1972) 257–273. [10] M.H. Kroll, J.D. Hellums, L.V. McIntire, A.I. Schafer, J.L. Moake, Platelets and shear stress, Blood 85 (1996) 1525–1541. [11] R.A. Adams, Sobolev Spaces, Academic Press, New York, 1975. [12] O. Ghattas, J.H. Bark, Optimal control of two- and three-dimensional Navier–Stokes flow, J. Computat. Phys. 136 (1997) 231–244. [13] M.D. Gunzburger, L.S. Hou, Treating inhomogeneous essential boundary conditions in finite element methods and the calculation of boundary stresses, SIAM J. Numerical Anal. 29 (1992) 390–424. [14] L.S. Hou, S.S. Ravindran, A penalized Neumann control approach for solving an optimal Dirichlet control problem for the for the Navier–Stokes equations, SIAM J. Control Optim. 36 (1998) 1795–1814. [15] V. Girault, P.A. Raviart, Finite Element Methods for the Navier–Stokes Equations, Springer, New York, 1986. [16] F. Brezzi, M. Fortin, Mixed and hybrid finite element methods, Computational Mathematics, vol. 15, Springer, New York, 1991. [17] T.E. Tezduyar, Y.J. Park, Discontinuity capturing finite element formulations for nonlinear convection-diffusion-reaction problems, Comput. Meth. Appl. Mech. Eng. 59 (1986) 307–325. [18] K.E. Jansen, S.S. Collis, C. Whiting, F. Shakib, A better consistency for low-order stabilized finite element methods, Comput. Meth. Appl. Mech. Eng. 174 (1999) 153–170.