A saddle point approach to an optimal boundary control problem for

In this paper we propose a saddle point approach to solve boundary control problems for the steady Navier-Stokes equations with mixed Dirichlet-Neumann boundary conditions, both in two and three dimensions. We provide a comprehensive theoretical framework to address (i) the well posedness analysis for the optimal control problem related to this system and (ii) the derivation of a system of first-order optimality conditions. We take advantage of a suitable treatment of boundary Dirichlet controls (and data) realized by means of Lagrange multipliers. In spite of the fact that this approach is rather common, a detailed analysis is still missing for mixed boundary conditions. We consider the minimization of quadratic cost (e.g., tracking-type or vorticity) functionals of the velocity. A descent method is then applied for numerical optimization, exploiting the Galerkin finite element method for the discretization of the state equations, the adjoint (Oseen) equations and the optimality equation. Numerical results are shown for simplified two-dimensional fluid flows in a tract of blood vessel where a bypass is inserted; to avoid to simulate the whole bypass configuration, we represent its action by a boundary velocity control.


Introduction
Optimal control problems (OCPs) for partial differential equations is a research area that features challenging theoretical questions for mathematical analysis, the need of devising efficient numerical algorithms, with a broad range of meaningful applications in engineering, like, e.g., aerodynamics [29,30] and biofluid dynamics [33,40].In this paper we address one of the most classical problems in this field, namely the boundary control for the Navier-Stokes equations.Much has been done for this problem, both on the mathematical analysis and on its numerical approximation; there are, however, some issues that have not yet been addressed -the case of the state problem with mixed Dirichlet and Neumann conditions is one of those.At the best of our knowledge, a comprehensive framework to address both (i) the well posedness analysis for the OCP related to this system and (ii) the derivation of a system of optimality conditions, is still missing.Setting this framework, and exploiting it in view of the numerical approximation of the OCP, is the scope of our paper.
More precisely, we study the mixed boundary value problem for the incompressible steady Navier-Stokes equations, both in two and three dimensions, where we prescribe homogeneous Neumann conditions (on the normal stress) on a portion of the boundary, homogeneous Dirichlet conditions on the velocity field on another boundary portion, and a velocity control on the remaining portion of the domain boundary * .For the sake of application, we consider the minimization of quadratic cost functionals involving either the square of the L 2 norm of the vorticity, or a velocity tracking-type term, on the whole computational domain.For the problem at hand, we introduce an ad-hoc saddle point formulation of the state problem which, besides velocity and pressure fields, requires a third variable representing the normal Cauchy stress and a Lagrange multiplier penalizing the Dirichlet data (this latter involving the control function, too).This saddle point formulation is more apt for numerical approximation than the one with a single field u that is obtained after projecting the Navier-Stokes equations onto a divergence free space.For the control problem, the saddle point formulation enables us to derive a system of first-order optimality conditions, and to carry out a rigorous mathematical analysis of the OCP.
Similar problems have been studied by several authors.What follows is an account of the most relevant works on the subject.In a classical paper by Gunzburger [21], the same control problem is addressed with however full Dirichlet boundary conditions (and a different cost functional), by a saddle-point formulation; a similar approach has been applied to viscous drag reduction in [19].In a series of other papers the problem is faced working on the subspace of divergence free velocity fields.In [12], the unsteady problem is studied, with Dirichlet conditions on the whole boundary; in [13], an exterior problem in two dimensions for the steady Navier-Stokes equations with boundary Neumann control is considered for the minimization of a drag functional; an inequality state constraint is also enforced on the velocity field.Distributed and Neumann controls are also considered, e.g., in [20].Boundary control problems for vorticity minimization have been analyzed in [32], in the case of steady Stokes equations, and in [31], considering instead unsteady Navier-Stokes equations; in both cases, however, a pure Dirichlet problem is considered.In [17] the analysis of existence of optimal boundary control for the Navier-Stokes equations with mixed boundary conditions is carried out in three dimensions, without deriving a system of optimality conditions.Regarding the case of unsteady Navier-Stokes equations, a boundary OCP similar to the one studied in our paper has been analyzed in [24], however dealing with a pure Dirichlet problem.The analysis and the approximation of OCPs involving distributed controls and in the case of unsteady Navier-Stokes flows has been considered in [22,23], again dealing with full Dirichlet boundary conditions.Mixed Dirichlet-Neumann boundary conditions have been instead considered in [4], however without carrying out the well-posedness analysis.We believe that our presentation is more comprehensive, and allows the reader to have a better picture of the whole mathematical analysis for both the state and the control problem and the derivation of a system of optimality conditions.
The structure of the paper is as follows.In Section 2 we introduce a ad-hoc saddle point formulation of the state problem with a Lagrange multiplier that represents the normal Cauchy stress on the whole boundary and carry out its well-posedness analysis in Section 3; some technical results related to the Stokes solver are reported in Section 4; In Section 5 we set up the ground for the analysis of the control problem; we then investigate the differentiability of the control-to-solution map in Section 6, and derive first order optimality conditions for the OCP in Section 7; In Section 8 we illustrate a possible way to approximate our control problem, by relying on a descent iterative method for the minimization of the cost functional, and using the Galerkin finite element method to discretize both the state and the adjoint problem; Numerical results shown in Section 9 concern the solution of an interior flow problem, involving a simplified two-dimensional tract of blood vessel where a bypass is inserted.The goal is to regularize the flow by minimizing a tracking-type quadratic functional, expressing the distributed observation of the (square of) the L 2 norm of the difference v − z d between the fluid velocity v and the target velocity z d , this latter being provided by the solution of the corresponding Stokes problem over the whole domain.Conclusions will then follow in Section 10.

Problem formulation
In this paper we consider a state system given by the steady Navier-Stokes equations, in velocity and pressure formulation, for steady, incompressible and viscous fluid flows.For this state system, we analyze a minimization problem for quadratic functionals of the velocity field through a velocity control operated on a portion of the inflow boundary.
Γ w denote the state space, endowed with the norm v X = ∇v 0 where Our goal is to minimize either the vorticity functional or a velocity tracking-type cost functional is a weak solution of the following system: where v and π denote the fluid velocity and pressure, respectively, ν is a given kynematic viscosity, and n is the unit outward normal direction.In the case of J2 , a flow matching problem is solved, aiming at driving the velocity field toward a target velocity profile z d ∈ X; this latter can be given, e.g., by the Stokes flow in the same domain, whenever the goal is to minimize recirculation and vortex generation for the sake of regularization of flows.If τ > 0, the second term in (2.1) or (2.2), called penalization term, besides helping in the wellposedness analysis, prevents using "too large" controls in the minimization of the cost functional.
We handle the Dirichlet boundary conditions on Γ 1 (involving the control function) by a Lagrange multiplier approach.Let e and b be the following bilinear forms and c be the trilinear form Using Hölder inequality, and the Sobolev inequality we have that Let us now define assuming that v in ∈ H 1/2 00 (Γ in ) d , we have that g ∈ H 1/2 00 (Γ 1 ) d , since u ∈ H 1 0 (Γ c ) d ⊂ T c with continuous injection, that is, for some positive constant c 0 , For the definition of these functional spaces, see, e.g., [7, Vol.III, Chap.VII].Observe that if v ∈ X, then v| Γ 1 ∈ H 1/2 00 (Γ 1 ) d .For ease of notation, hereon we denote by T 1 = H 1/2 00 (Γ 1 ) d and by •, • Γ 1 the duality pairing between T 1 and its dual T * 1 .In the same way, we will denote by T c = H 1/2 00 (Γ c ) d , and by •, • Γ c the duality pairing between T c and its dual The weak formulation of (2.3) that we consider reads as follows: find (2.7) Note that a formal integration by parts in the first equation gives for the boundary stress πn − ν ∂v ∂n , and 0 = πn − ν ∂v ∂n on Γ out .
Thus, the first equation in (2.7) enforces in a natural way the homogeneous Neumann condition on Γ out .By grouping together the last two equations in (2.7), we obtain the following equivalent formulation: Upon defining the bilinear form b : we reformulate problem (2.9) as: (2.11)

Analysis of the state problem
Despite their ubiquitous use when modeling viscous incompressible flows, Navier-Stokes equations endowed with mixed boundary conditions have been seldom analyzed from a theoretical standpoint.If homogeneous Neumann conditions are set on a boundary portion, energy inequalities or a priori estimates of a weak solution to the Navier-Stokes system cannot be obtained because there is no control over the energy flux on this boundary portion.Some results can be found, e.g., in [41] regarding the existence and uniqueness of solution in weighted L p Sobolev spaces and Hölder spaces for problems set over polyhedral domains; see also, e.g., [3,35,36] for the time-dependent case and [2] for a wellposedness analysis in divergence-free spaces.Additional conditions that estimate backward flows over the Neumann boundary have been considered, e.g., in [34].
In this paper we derive a well-posedness result for the steady Navier-Stokes system in presence of mixed Dirichlet-Neumann boundary conditions including the pressure field among the state variables, in non-weighted Sobolev spaces.We show that if g T 1 is small enough, then (2.10) is well posed.We first need to introduce the following continuous operators, related to the bilinear form b: We denote by N B the kernel of B, by N ⊥ B its orthogonal (in X) and by N 0 B the polar of N B , that is The so called inf-sup condition holds for the bilinear form b: Lemma 3.1.There exists a constant β > 0 such that Proof.Given ( φ, λ) ∈ Q, according to Riesz Theorem there exists a unique λ * ∈ T 1 such that . Such ṽ can be obtained as follows.Let r ∈ X be the lifting of Let us now set ṽ = w + r, with w ∈ H 1 0 (Ω) d to be determined.We have It follows (see, e.g., [15]) that there exists w ∈ H 1 0 (Ω) d with divw = φ − divr so that ṽ = w + r meets the required assumptions.Thus we can write b(ṽ, ( φ, λ)) , by taking the supremum over ṽ and the infimum over ( φ, λ).The following lemma is well known (see [15,Chap. 1,Thm. 4.1]).Lemma 3.2.Condition (3.1) is equivalent to each one of the following properties: We are ready to prove the following result.
To solve problem (2.10), we first determine v 0 by solving where G ∈ Q * is defined in (2.11).Then we determine z ∈ N B from the equation where and that F, We divide the proof in four steps.
Step 1. (3.10) Step 2. Now we reduce (3.7) to a fixed point equation.Let z ∈ N B be given and consider the linear For each z, the bilinear form ẽ (z; with η as in (3.5), we have and, for m integer, define By the Lax-Milgram Theorem, for each z ∈ S 4δ u , we can define a solution map operator Thanks to the continuity of ẽ, we can rewrite equation (3.15) as where A (z) ∈ L(N B , N * B ). Clearly T (z) = A (z) −1 and therefore we have the identity L is a strict contraction and therefore has a unique fixed point in S δ . Moreover, Thus, T (z) F 0 : S 4δ u → S 4δ u is a strict contraction and has a unique fixed point Step 4. It remains to solve equation (3.8).Observe that equation (3.7) means that the operator E (v) ∈ X * given by where v = z * + v 0 , belongs to N 0 B .Therefore, by Lemma 3.2, (ii), there exists a unique (π, t) and the proof is complete with C 1 = 5 β−1 and C 2 depending only on the constants ν, β, C S .
Remark 3.5.Theorem 3.3 holds with a right hand side f ∈ X * with small enough f X * norm.Then the unique solution (v,π, t) satisfies the estimates where C 1 and C 2 only depend on the constants ν, β, C S .
Under the hypotheses of Theorem 3.3, stressing the dependence on u of the solution to problem (2.10), we define the solution map S : u → (v (u) , π (u) , t (u)) .We will need the following result, concerning the map u → v (u), since only this component is involved in the cost functional.Proposition 3.6.Under the hypotheses of Theorem 3.3, possibly with a smaller η depending only on ν, β, C S , there holds with C L depending only on ν, β, C S .
From (3.10), since the equation for v 0 j is linear and The equation for z = z 2 − z 1 reads, after routine calculations, Similarly, using (3.21) and v 0 j X ≤ 5 β−1 η,

Collecting all the above inequalities, and noting that
then (3.19) follows with C L depending only on ν, β, C S .

The Stokes solver
The same technique used in the proof of Theorem 3.3 provides existence and uniqueness of the weak solution (v,π, t) ∈ V of the Stokes problem for any where C only depends on the constants ν, C S .Let us introduce the Stokes solver is the unique solution of (4.1).By (4.2), R is a (linear) continuous isomorphism between V * and V.
Moreover, the following result holds.

Existence of an optimal control
We are now ready for the well-posedness analysis of the optimal control problem Here J is given by either (2.1) or (2.2) and , where η denotes the constant defined in (3.5).Note that the definition of U ad depends on the positive constant c 0 .If we assume that there exists an optimal control û ∈ U ad and a corresponding optimal state v ∈ X, (π, t) ∈ Q.
and recall that X is compactly embedded into Finally, thanks to the weak convergence of {v k } in X, we have Passing to the limit for k → ∞ into (5.3),we deduce that (v, π, t, u) is a solution of (2.7).(4) J is continuous and convex and therefore also sequentially weakly lower semicontinuous.

Differentiability of the solution map
In order to derive first order optimality conditions, we need to show that the solution map S : u → (v (u) , π (u) , t (u)) is Fréchet differentiable in U ad .Precisely, we show the following result.
The solution map S is Fréchetdifferentiable at ū and its derivative S ( ū) is given by where, setting v = v ( ū), (w, q, s) ∈ V is the solution to (6.1) Remark 6.2.Since ū X < η, the bilinear form a (w, ϕ) = e(w, ϕ) is X−coercive hence the linear problem (6.1) is well posed.
Remark 6.3.From the previous proof, it is clear that S is differentiable also with respect to the T cnorm.

First order optimality conditions
We are now ready to write the first order optimality conditions for our OCP (5.1).Precisely, we can show the following result.Theorem 7.1.Let û be an optimal control and (v, π, t) the corresponding optimal state.There exists a unique multiplier (ẑ, q, ŝ) ∈ V such that û and (ẑ, q, ŝ) satisfy the following optimality system: • variational inequality: for every k ∈ U ad , where J(u) = J(v(u), u) is the reduced cost functional.
Proof.Define the nonlinear operator N : The nonlinear component of N is a continuous bilinear operator and therefore N is Fréchet differentiable ‡ , with and N u (v, π, t, u) k = (0, 0, −g 0 (k)) (7.4) where Consider now the Stokes solver R : V * → V defined in Section 4. In terms of R and N, the state problem can be written as In Lemma 7.3 below we prove that: (a) the adjoint operator N (v,π,t) (v, π, t, u) * : : V → V * is given by where n is the exterior unit normal to Γ c ; (b) the adjoint operator N u (v, π, t, u) * : V → T * c is given by N u (v, π, t, u) * : (z,q, s) = −s. (7.7) By a well known theorem (see, e.g., [27,Thm. 1.45]) the existence of a unique multiplier (ẑ, q, ŝ) ∈ V follows from the properties (i) , (ii) , (iii) listed below.
(i) J and G are Fréchet differentiable, with (ii) From Theorem 6.1 we know that, for every u ∈ U ad , the state equation (7.5) defines a unique solution map S (u) which is Fréchet-differentiable.
in weak form reads: (7.9) Since ū ∈ U ad , the bilinear form is X−coercive, therefore problem (7.9) is well posed and G (v,q,s) v, π, t, ū is a continuous isomorphism between V and V * .Introducing the Lagrangian functional where (z,q, s) is a multiplier, the first order optimality conditions are obtained by differentiating L.
Precisely, the adjoint problem reads ( L or, equivalently, upon introducing the adjoint operator for every (ψ, η, r) ∈ V.
As in Remark 6.2, since û ∈ U ad , the bilinear form (adjoint of the bilinear form (6.2)) is X−coercive and the linear problem (7.1) is well posed.To write the adjoint system in strong form, note that Moreover, in the case of the vorticity functional J1 , Therefore, the adjoint system in strong form reads, setting the rotation tensor ∇v − (∇v We can also deduce that, regarding the multiplier, In the case of the velocity tracking-type cost functional J2 , the strong form of the adjoint problem is straightforward to obtain, and reads as follows: In this case, we find instead that ŝ = qn − ν ∂ẑ ∂n on Γ 1 .(7.17) Remark 7.2.From (7.2) we deduce that the derivative of the reduced functional J (u) = J (v (u) , u) at u = û ∈ U ad is given by Then, introducing the embedding operator B : The variational inequality (7.2) can be written in the form If Bŝ 0 and τ = 0, (7.19) gives If Bŝ 0 and τ > 0, we have û = − ŵ/τ if ŵ/τ ∈ U ad , otherwise û is given by (7.20).
Lemma 7.3.The following expressions hold: where n is the exterior unit normal to Γ c .

Numerical approximation
For the numerical solution of the OCP (5.1) we use an iterative method for the minimization of the cost functional, and we exploit the Galerkin finite element method for the numerical discretization of the state equation, the adjoint equation and the optimality equation.In this way we can take advantage of the system of first-order (necessary) optimality conditions we have derived in the previous sections.Several alternatives have been proposed in literature; see, e.g., [14,25,28,43] for the case of sequential quadratic programming methods applied to boundary control of steady Navier-Stokes flows, [9] for an augmented Lagrangian method, and [5] for inexact Newton methods.A review of numerical strategies for OCPs can be found, e.g., in [1,18,26].

Numerical discretization of state, adjoint and optimality equations
The state system consists of the steady Navier-Stokes equations, while the adjoint system is given by a steady Oseen problem.Following the saddle-point approach described so far, for both these problems Dirichlet boundary conditions are enforced in weak form, by means of the Lagrange multipliers t and s, respectively.These two latter variables are indeed related with the boundary stress over the Dirichlet boundary, see, e.g., equations (2.8) and (7.15)-(7.17).As a matter of fact, for the numerical approximation, we can rewrite the state and the adjoint equations in a conventional manner, by introducing appropriate functional spaces and suitable lifting terms to take into account non-homogeneous Dirichlet conditions.This option allows us to avoid the explicit computation of the Lagrange multipliers t and s.
For the numerical approximation of both these problems, let us introduce two subspaces X h ⊂ X and Q h ⊂ L 2 (Ω), of dimension N V , N Q < +∞, respectively, and let us denote by v h ∈ X h and p h ∈ Q h the FE approximations for the velocity and the pressure fields [44].Similarly, z h ∈ X h and q h ∈ Q h will denote the FE approximations for the adjoint velocity and pressure fields, respectively.We also denote by u h ∈ U h ∩ U ad the FE approximation of the control variable, where U h ⊂ U is a subspace of the control space U = H 1 0 (Γ c ) d of dimension N U < +∞.The Galerkin-FE approximation of the state system (2.7) thus reads as follows: we seek for where All the forms are continuous over the discrete spaces X h and Q h .Conversely, the approximation stability is ensured by imposing that the coercivity and inf-sup conditions are still valid at the discrete level.In particular, the coercivity of e over X h is inherited from that over X; on the other hand, we require that b is inf-sup stable over X h × Q h , so that the following discrete Brezzi inf-sup condition, holds [6]: ∃ βh > 0 : This last property is ensured e.g., by choosing X h × Q h as the space of Taylor-Hood P 2 − P 1 finite elements for the velocity and the pressure, respectively (in both the state and the adjoint problem); however, this choice is not restrictive -the whole construction keeps holding for other spaces combinations as well [16].
We can now derive the matrix formulation corresponding to the Galerkin-FE approximation (8.1).Let us denote by {φ v i } N v i=1 and {φ p i } N p i=1 the Lagrangian basis of the FE spaces X h and Q h , respectively, so that we can express the FE velocity and pressure as We remark that the solution to (8.1) is vanishing on the whole Dirichlet boundary, so that the corresponding velocity approximation fulfilling the boundary conditions is given by v h + v D h .Thus, by denoting as v h ∈ R N v and p h ∈ R N p the vectors of the degrees of freedom appearing in (8.3), problem (8.1) can be rewritten as: where, for 1 ≤ i, j ≤ N v and 1 We solve the nonlinear saddle-point problem (8.4) by means of a fixed-point method (also referred to as Picard iteration, see, e.g., [11,Sect. 8.2] it has a larger ball of convergence than Newton's method (see e.g., [11], Chapter 7.2 and references therein).Moreover, if a small data condition (like, e.g., (3.3) for the case at hand) is satisfied, the fixed-point method is globally convergent.Thus, starting from an initial guess (v (0)  h , p (0) h ), for k ≥ 1 we solve ) given a small tolerance ε NS tol > 0. As initial guess, we take the Stokes solution of (8.4).Each Oseen system (8.6) is solved by means of a sparse LU factorization.
Regarding instead the numerical approximation of the adjoint variables, the adjoint system (7.11) is linear with respect to both the adjoint velocity and pressure.Moreover, it depends on the state variables only through the term c (v, ψ, ẑ) + c (ψ, v, ẑ) and the right-hand side of (7.11); hence, the FE arrays corresponding to these two terms must be reassembled at each step during the optimization procedure.Given the approximation v h = v h (u h ) of the (state) velocity and u h of the control function, the Galerkin-FE approximation of the adjoint system (7.11)thus reads as follows: we seek for (z where

Numerical optimization
We now describe the numerical optimization method exploited for the solution of our OCP.We rely on a gradient-based algorithm exploiting the gradient ∇J (u h ) to iteratively update the control until a suitable convergence criterion is fulfilled; a similar approach can be found, e.g., in [8].Notable instances are descent methods, such as the gradient, (non-linear) conjugate gradient, quasi-Newton or Newton methods.In the simplest case of a gradient (or steepest descent) method, starting from an initial guess u (0) , we iteratively generate a sequence where λ (k) > 0 is a step size, until, e.g., a suitable stopping criterion is fulfilled, and P U ad denotes the projection onto U h ∩ U ad , to be evaluated similarly to what we have done in (7.20).
As stopping criterion, we might require that either for a given tolerance ε > 0. At each step of the gradient method, we solve the state Navier-Stokes system (8.1), from which we obtain v (k) , and the value of the cost functional J(v (k) , u (k) ).Then, we solve the adjoint system (8.7)(v (k) and u (k) being given) and obtain (z (k) , q (k) ).From these variable, we determine the boundary stress and approximate (by the Galerkin finite element method) the solution of the following problem, We highlight that at each step of the steepest descent method, we take as initial guess of the Newton method used to solve the state Navier-Stokes system the state velocity computed at the previous step of the optimization procedure, instead of computing it by solving a Stokes system, for the sake of computational efficiency.

Numerical results
In this section we present the numerical results related with an OCP for the optimal control of blood flows -described through a steady Navier-Stokes model -through a bypass graft.This latter usually provides blood flow through an alternative bridging path in order to overcome critically occluded arteries; one of the most dangerous cases is related to coronary arteries, which supply the oxygen-rich blood perfusion to the heart muscle.The design of the connection between the graft vessel and the host arterial vessel is a critical factor in avoiding post-operative recurrence of the occlusion, since fluid dynamic phenomena such as recirculation, oscillating or untypical high/low shear rates, and stagnation areas, can cause the growth of another stenosis downstream the arterial-graft connection.Hence, a typical design of a bypass graft aims at minimizing some cost functionals related to haemodynamic quantities [10,42].However, a rigorous model for blood circulation should take into account (i) the flow unsteadiness, (ii) the arterial wall deformability, and possibly (iii) complex rheological model to characterize the aggregated nature of blood.Optimal control (and optimal design) problem require the repeated simulation of these flow equations (and the evaluation of the cost functional to be minimized); to reduce the computational burden, we adopt steady incompressible Navier-Stokes equations for laminar Newtonian flows.
We focus on the minimization of the tracking-type functional (2.2) in order to drive the blood velocity towards a specified velocity target state z d , featuring a regular pattern; this latter is given by the solution of the Stokes system in the same domain.Other possible cost functionals for the problem at hand have been employed, e.g., in [37,38]; the effect of uncertainty (e.g., affecting either the residual flow across Γ in , or the geometrical configuration of the bypass) has been explored in [38,39,45].
Here a two-dimensional geometrical configuration is considered, although the whole theoretical analysis and the computational pipeline are perfectly valid for three-dimensional problems as well.
For the case at hand, we consider a fluid flow in a tract of blood vessel where a bypass is inserted.As a matter of fact, the whole bypass graft is not simulated; its action is represented via a velocity control u acting on the boundary Γ c ⊂ ∂Ω, the interface where the final portion of the graft and the host vessel meet (see Figure 1).We consider two different scenarios: (i) a total occlusion of the host artery (setting v in = 0 on Γ in , and (ii) the presence of a residual blood flow (v in 0) on that boundary.We choose ρ = 1, v = 25cm s −1 , d = |Γ c | ≈ 0.75cm from which, defining the Reynolds number Re = ρ v d/µ, which is set equal to Re = 200, the dynamic viscosity coefficient µ can be deduced.We first compute the solution to the Stokes system, to determine the target velocity profile z d ; we set the relaxation parameter λ (k) = τ = 0.1 and the regularization parameter τ = 0.05.We consider the stopping criterion based on the difference of two successive iterates (8.10), for which we select ε = 0.01; the tolerance for the stopping criterion of the Newton method is instead given by ε NS tol = 10 −5 .For our computations we use a mesh with 27, 233 triangular elements and 14, 048 vertices.
We first consider the case of total occlusion of the host arterial vessel, that is, v in = 0.In Figure 2 we report the velocity (both its magnitude and the streamlines) and pressure fields corresponding to the initial control function u (0)  h and the optimal control ûh obtained when the steepest descent algorithm stops, after N it = 157 iterations.In this case, the cost functional is reduced of about 37%, decreasing from J(u (0) h ) = 288.241 to J( ûh ) = 181.445.A reduction in the magnitude of the adjoint velocity field, which can be remarked by comparing the plots displayed in Figure 3, also shows that the derivative of the cost functional evaluated at u h = ûh is smaller than the corresponding quantity calculated when

Mathematics in Engineering
Volume 1, Issue 2, 252-280.u h = u (0) h .The magnitude of the vorticity is reported in Figure 4: a global decrease of this quantity is indeed obtained also by minimizing the tracking-type cost functional, as it results by comparing the vorticity magnitude for the initial and the optimal control.actually reduces the recirculation around the graft-artery junction, (compare Figure 5 and Figure 2).Also in this case a remarkable reduction of the value of the cost functional is obtained: the steepest descent algorithm stops after N it = 151 iterations, and the cost functional is reduced of about 37%, decreasing from J(u (0) h ) = 289.796 to J( ûh ) = 183.565.Similar considerations can be obtained by comparing the adjoint velocity and the vorticity fields calculated for the initial and the optimal controls (Figures 6 and 7).
For a concrete example we denote by Ω ⊂ R d , d = 2, 3 the domain in Figure 1 below, occupied by the fluid and by Γ = ∂Ω its boundary.The inflow boundary Γ 1 = Γ c ∪ Γ in is composed by Γ c , the control boundary, and by Γ in , where a Dirichlet data is assigned; Γ out and Γ w are the outflow and the wall boundary portions, where zero stress and zero velocity are assigned, respectively.

Figure 1 .
Figure 1.The domain and its boundary portions.

Figure 2 .
Figure 2. State velocity magnitude (top) and pressure with velocity streamlines (bottom) for the initial control (left) and the optimal control (right).Case v in = 0.

Figure 3 .
Figure 3. Adjoint velocity magnitude for the initial (left) and the optimal (right) control.Case v in = 0.

Figure 4 ., 0 T
Figure 4. Vorticity magnitude for the initial (left) and the optimal (right) control.Case v in = 0. Similar results are obtained when a residual flow across Γ in is present; in this case, we set v in = a exp − (y − ȳ) 2 2b 2 , 0 T

Figure 5 .
Figure 5. State velocity magnitude (top) and pressure with velocity streamlines (bottom) for the initial control (left) and the optimal control (right).Case v in 0.

Figure 6 .
Figure 6.Adjoint velocity magnitude for the initial (left) and the optimal (right) control.Case v in 0.

Figure 7 .
Figure 7. Vorticity magnitude for the initial (left) and the optimal (right) control.Case v in 0.
45. Sankaran S, Marsden A, (2010) The impact of uncertainty on shape optimization of idealized bypass graft models in unsteady flow.Phys Fluids 22: 121902.c 2019 the Author(s), licensee AIMS Press.This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)