Error Bounds for Discontinuous Finite Volume Discretisations of Brinkman Optimal Control Problems

We introduce a discontinuous finite volume method for the approximation of distributed optimal control problems governed by the Brinkman equations, where a force field is sought such that it produces a desired velocity profile. The discretisation of state and co-state variables follows a lowest-order scheme, whereas three different approaches are used for the control representation: a variational discretisation, and approximation through piecewise constant and piecewise linear elements. We employ the optimise-then-discretise approach, resulting in a non-symmetric discrete formulation. A priori error estimates for velocity, pressure, and control in natural norms are derived, and a set of numerical examples is presented to illustrate the performance of the method and to confirm the predicted accuracy of the generated approximations under various scenarios.


Introduction
Fluid control problems are highly important in diverse fields of science and engineering. For example, they are encountered in the minimisation of drag, in the design of devices serving to increase mixing properties, in the reduction of turbulent kinetic energy, and in several other applications. Some of the earliest references to theoretical aspects of these control problems can be found in the classical works [1,35]. The literature that relates to their numerical approximation is quite abundant, especially if associated to finite element methods (see e.g. [9,26,45,47,48] and the references therein). Works focusing on the approximation of control problems subject to Stokes and Navier-Stokes flows typically employ conforming discretisations for state, co-state and control variables. It has been found that the convergence rate of the control approximation is of O(h) and O(h 3 2 ) for piecewise constant and piecewise linear discretisations, respectively. On the other hand, using the so-called variational discretisation approach (cf. [29], in which the control set is not discretised explicitly but recovered by a projection), it is possible to improve this convergence rate to O(h 2 ). Alternatively, a similar behaviour is observed if one uses graded meshes instead of uniform partitions [44], or if using piecewise constant control discretisations when state and adjoint variables are approximated with Lagrangian finite elements [42].
The present paper focuses on finite volume element (FVE) approximations (methods where one introduces a dual mesh and reformulates a pure finite volume scheme in the form of a Petrov-Galerkin scheme). A priori error estimates for FVE schemes applied to linear elliptic and parabolic optimal control problems have been established in [38,39]. These methods are based on the optimise-then-discretise approach, which we will adopt herein. In this context, we recall that the order in which the optimisation and discretisation steps are performed, results in different discrete adjoint equations and the solutions may not coincide (see the review [9] and the references therein). We will concentrate the analysis on a particular class of FVE schemes: a hybrid strategy called discontinuous finite volume (DFV) method, where discontinuous piecewise linear functions conform the trial space, and piecewise constant test functions are used in a FV fashion. The application of these schemes in the approximation of Stokes and related fluid problems can be found in e.g. [10][11][12]24,32,33,54].
Discontinuous approximations will be generally preferable to guarantee preservation of physically relevant properties. They would also be appropriate when the model exhibits rough coefficients and where sharp solutions are expected. Permeability fields possess this behaviour in many applicative scenarios, and DFV schemes would be of special interest. Other advantages of DFV formulations include flexibility for choosing accurate numerical fluxes, smaller dual control volumes, and suitability for error analysis in the L 2 -norm. In the formulation advanced herein the momentum equation is tested against vector functions spanned by a basis associated to a dual grid, and the mass conservation equation is tested against piecewise constants defined on the primal mesh. Integration by parts on each dual element yields a classical finite volume scheme defined in terms of fluxes across the boundaries of dual elements. Then, some particular features of a given lumping map connecting discrete functions associated with the primal and dual meshes allow us to rewrite the formulation completely in terms of volume integrals involving primal elements, except for the zeroth-order term, the right-hand sides of the state and costate equations, and all jump terms that appear in the off-diagonal bilinear forms.
For non-viscous flow in porous media written in primal (pressure) formulation, the permeability tensor manifests itself as an anisotropic diffusion, and some methods are available for their successful discretisation. These include cell and vertex-centred schemes [23,55], discrete duality finite volumes [16], high-order gradient reconstruction finite volume methods [40], mimetic schemes [18], anisotropic FVE methods [36], virtual finite volumes [22], or discontinuous immerse FVE schemes [37]. Even if the treatment is simpler in our case, as the inverse permeability is assumed isotropic, and it appears only in the drag term; our intention is not to perform a thorough comparison against these techniques, but rather to regard our contribution as a natural extension of the optimal control problems solved using specifically DFV methods (and having so far being constructed for systems governed by linear and semilinear elliptic, semilinear parabolic, and hyperbolic equations [49][50][51][52]) to the case of velocity control for the Brinkman equations. We emphasise once more that the discontinuous character of permeability fields represents a clear motivation for employing DFV methods. On the other hand, for the approximation of the control variable, we will discuss three alternatives: a variational discretisation approach, element-wise constant and element-wise linear discretisation.
The paper is structured in the following manner. The remainder of this section includes some standard notations, statement of the governing problem along with its weak formulation, and the corresponding optimality condition in continuous form. Next, in Sect. 2 we formulate the DFV scheme of the considered optimal control problem. Section 3 focuses on the development of a priori error estimates for different types of control discretisations. Finally, in Sect. 4 we summarise the solution algorithm and illustrate our theoretical error bounds and performance of the method by a set of numerical experiments.
Let Ω ⊂ R d , d = 2, 3, be a bounded convex polygonal domain with boundary ∂Ω. The outward unit normal vector to Ω is denoted by n. Standard terminology will be employed for Sobolev spaces: The corresponding norms will be denoted by · 1,Ω . We also consider the space of integrable functions with zero mean: L 2 0 (Ω) = {q ∈ L 2 (Ω) : Ω q dx = 0}, and will write L 2 (Ω) = L 2 (Ω) d . By div we will denote the usual divergence operator div applied row-wise to a tensor, I will denote the d × d identity matrix, and 0 will be used as a generic null vector. The Optimal Control Problem Let us consider the following distributed optimal control problem min u∈U ad governed by the Brinkman equations where U ad is the set of feasible controls, defined for −∞ ≤ a j < b j ≤ ∞, j = 1, . . . , d by This model describes the motion of an incompressible viscous fluid within an array of porous particles, and according to the flow regime characterised by the ratio between permeability and viscosity, it can represent both the Darcy and Stokes limits. Here y denotes the fluid velocity, p is the pressure field, u is the control variable, and λ > 0 is a given Tikhonov regularisation (or control cost) parameter. The quantity νε( y)− pI is the Cauchy (true stress) tensor, where ε( y) = 1 2 (∇ y + ∇ y T ) is the infinitesimal rate of strain, ν(x) is the dynamic viscosity of the fluid, and K(x) stands for the permeability tensor of the medium divided by the viscosity. The forthcoming analysis requires that this matrix is isotropic. Here the desired velocity y d and the applied body force f are known data with assumed regularity L 2 (Ω) or H 1 (Ω), depending on the specific case. One seeks to identify an additional source u giving rise to a velocity y in order to match a target velocity y d . We stress that by proceeding analogously to the proof of [53,Theorem 2.37] it can be shown that the optimal value of the control u is in H 1 (Ω) under the assumption that y d ∈ L 2 (Ω). A similar observation has been made in [14], and used in the derivation of error estimates.
We assume that K is symmetric, uniformly bounded and positive definite, i.e., there exist two positive constants k 1 and k 2 such that We also assume that the variable viscosity satisfies The weak formulation associated to the state equations (1.2)-(1.4) is given by: where the bilinear forms a(·, ·) : for all y, v ∈ H 1 0 (Ω) and q ∈ L 2 0 (Ω). Above (·, ·) 0,Ω stands for the scalar product in L 2 (Ω) and · 0,Ω denotes the associated norm. The bilinear form b(·, ·) relating the functional spaces for velocity and pressure satisfies the following Babuška-Brezzi condition (see [46], for example): there exists ξ > 0 such that and, together with the ellipticity of a(·, ·)+c(·, ·), it implies the unique solvability of problem (1.6). The optimal control problem (1.1)-(1.4) under consideration is strictly convex. Hence, it admits a unique optimal solution, and the first order necessary conditions are also sufficient for optimality (for details on well-posedness and first order optimality we refer to [35]). The optimality condition can be formulated as J (u)(ũ − u) ≥ 0 ∀ũ ∈ U ad , and also rewritten as Fig. 1 Left: sketch of a single primal element T in T h , and sub-elements T * i belonging to the dual partition T * h . Right: its three-dimensional counterpart, showing a tetrahedron T decomposed into four sub-tetrahedra where w is the velocity associated with the adjoint equation (1.10) In turn, the variational inequality (1.7) can be equivalently recast in component-wise manner where the operator P denotes a projection defined for a generic scalar function g as P [a,b] (g(x)) = max(a, min(b, g(x))), a.e. in Ω.
It is not difficult to see that this projection satisfies the following regularity property (see also [43]) In addition to T h (from now on, referred to as primal mesh), we introduce a dual partition in the following way. Each element T ∈ T h is split into three sub triangles (or four subtetrahedra if d = 3) T * i , i = 1, . . . , d + 1, by connecting the barycentre of the element to its corner nodes (see a schematic for d = 2 and d = 3 in Fig. 1). The set of all these elements generated by barycentric subdivison will be denoted by T * h and will be called the dual partition of T h . Let e be an interior face shared by two elements T 1 and T 2 in T h . By n 1 and n 2 we will denote unit normal vectors on e pointing outwards T 1 and T 2 , respectively. The average {{·}} e and jump [[·]] e operators defined on e for a generic scalar or vectorial field v are: Note that jump and averages are defined so that they preserve the dimension of the argument.
We denote by P m (T ) the space of polynomials of degree less or equal than m, defined on T ∈ T h , and P m (T ) will denote its vectorial counterpart. A finite dimensional trial space (used for the state and co-state velocity approximation) associated with the primal partition T h is while the finite dimensional test space for velocity (and corresponding to the dual mesh T * h ) is Moreover, the discrete space for state and co-state pressure approximation is defined as and we define a space with higher regularity These spaces, associated with the two different meshes, are connected through the transfer operator γ : V(h) → V * h , characterised in the following manner: Some useful properties of this map are as follows.
Lemma 1 Let γ be a transfer operator defined as in (2.1). Then i) γ is self-adjoint with respect to the L 2 -inner product, i.e.
ii) The operator γ is stable with respect to the norm · 0,Ω , that is ,Ω , then |||·||| 0,h and · 0,Ω are equivalent, with equivalence constants being independent of h. iii) There exist constants C 0 and C 1 independent of h such that Proof The proof of (i) and (ii) are given in [7] for the scalar version of the inter-mesh operator (2.1), and the same arguments can be used to extend their validity to the vector case. Property (iii) follows analogously to the proof of [20, Lemmas 2.3 and 2.6]. For a proof of (iv) we refer to [31].
Let us stress that throughout the paper, the symbol C will represent a generic positive constant independent of meshsize h, that can take different values at different instances.

Discrete Formulation for the State and Adjoint Equations
Let v h ∈ V h . We proceed to test (1.2) and (1.3) against γ v h ∈ V * h and φ h ∈ Q h , respectively, and after integrating by parts the momentum equation on each dual element, and the mass conservation equation on each primal element, we end up with the following scheme: find [11]): where, A d+2 = A 1 . With the help of Fig. 1 we see that, in 2D, for a fixed j the integrals A j+1 B A j are considered over a path of two segments, whereas in 3D they are taken over a triangular facet. In any case, they contribute to construct the normal fluxes on the interior faces of each dual sub-element T * j , and so the symbol n also denotes the outer normal on that sub-triangle or sub-tetrahedron. In addition, the constants α d and δ = (d − 1) −1 are parameters independent of h commonly used in interior penalty methods.
Proceeding analogously, we can write down a DFV formulation for the adjoint equation (1.8)-(1.10) as follows: For the sake of our forthcoming analysis, we introduce the following discrete norms in V(h), which are naturally associated with the bilinear form c h (·, ·): and we note that they are equivalent in V h . Moreover, we also have the following discrete Poincaré-Friedrichs inequality (see [11, pp. 457 and we can use Cauchy-Schwarz inequality and the definition of γ to readily obtain Proceeding analogously to [56, Lemma 6], we can establish the coercivity of the bilinear form A h (·, ·), stated in the following result.

Lemma 2 Let us assume that
Proof Let B = K −1 and consider its average tensor, defined on each primal element by Then, in view of Cauchy-Schwarz inequality together with properties (2.9) and (2.2), we can assert that and therefore, relations (2.10) and (2.3) lead to

Lemma 3
The bilinear forms defined in (2.6) possess the following properties: and it satisfies the bound

ii) For the non-symmetric bilinear form c h (·, ·) it holds that
where for (2.12), α d > 0 is assumed sufficiently large. iii) The choice of approximation spaces V h and Q h yields the condition Proof For i) it suffices to apply the definition of A h (·, ·), together with relation (1.5), and the norm equivalence between |||·||| 0,h and · 0,Ω . Results in ii) have been established in [33] and [7], whereas proofs for iii)-iv) can be found in [54].

Discretisation of the Control Variable
Let U h ⊆ L 2 (Ω) denote the discrete control space, and let us introduce the discrete admissible space for the control field as U h,ad = U h ∩U ad . Three approaches are outlined in what follows.

Variational Discretisation
In the so-called variational approach (cf. [29]), control variables are not discretised explicitly, that is, one simply takes U h = L 2 (Ω) and in this case the discrete and continuous admissible spaces U h,ad and U ad coincide. Consequently, the control variable does not necessarily lie in a finite element space associated to T h , and typically one requires a nonstandard implementation and more involved stopping criteria for the algorithms of control computation. Discretisation errors using this method will be addressed in Sect. 3.2.

Piecewise Linear Control Discretisation
Here we approximate the control variable with the similar elements as those employed for state and co-state velocity. That is, It is worthy to note that the state velocity space V h coincides with the control space in the case of homogeneous Neumann boundary conditions, whereas for Dirichlet boundary data,

Piecewise Constant Discretisation
In this case, the discrete control space is defined as The convergence properties associated with the above two approaches will be derived in Sect. 3.3, but already at this point we can apply Lemma 3 along with the Babuška-Brezzi theory for saddle point problems to ensure the unique solvability of (2.4)-(2.5), for a fixed control u h .

Convergence Analysis
In this section we provide a priori error estimates for DFV approximations of the state and adjoint equations, and for the three control discretisation approaches outlined in Sect. 2.3.

Preliminaries
For a given control u and f , let the pair ( y h (u), p h (u)) be the solution of the following problem Similarly, for a given state velocity y, let (w h ( y), r h ( y)) be the solution of We then proceed to decompose total errors in the following manner:

Lemma 4 There exists a positive constant C independent of h such that the following estimates hold
Adding (3.9) and (3.10) after In turn, using the coercivity of A h (·, ·) and c h (·, ·) in combination with (2.2) and (2.7), we obtain which readily yields the bound On the other hand, applying the inf-sup condition (2.14), using (3.9), the boundedness of A h (·, ·) and c h (·, ·), along with (3.11), we realise that (3.14) Proof We proceed analogously to the proof of [24, Theorem 3.1] and directly apply Lemma 3 to readily derive the following estimates: Next, the derivation of L 2 -estimates for y − y h (u) and w − w h ( y) follows an Aubin-Nitsche duality argument. Let us consider the dual problem: which is uniquely solvable, and moreover the following H 2 (Ω) × H 1 (Ω)-regularity is satisfied: where the auxiliary bilinear forms adopt the following expressions Since z I ∈ V h is a continuous interpolant of z, we note that the pair y− y h (u), p− p h (u) will be a solution of the following problem Using the definition of c h (·, ·) and C h (·, ·) we can assert that where the inner product over the primal mesh is understood as the sum of the inner products over each element in T h . On subtracting Eq. (3.19) from the sum of Eqs. (3.18) and (3.20), and using (3.21), it follows that Notice that the estimation of R 1 results as a combination of the boundedness of K −1 , assumption (1.5), the bounds (3.17), the self-adjointness and approximation properties of γ stated in (3.16), and Cauchy-Schwarz inequality. This gives where the last inequality follows from (3.13). For the second term we employ the definition of c h (·, ·), and relations (3.17), (3.16) to verify that Bounds for the remaining terms can be obtained following the proof of [33,Theorem 3.4] and [24, Theorem 3.2], as follows Combining the five estimates above with (3.22), we straightforwardly obtain and very much in the same way, one arrives at w − w h ( y) 0,Ω ≤ Ch 2 w 2,Ω + r 1,Ω + y 1,Ω + y d 1,Ω .
Now, for a given control u, let (w h (u), r h (u)) be the solution of and notice that similar arguments as those appearing in the proof of Lemma 5 and in the derivation of the estimate y − y h (u) 0,Ω ≤ Ch 2 , will readily lead to The following result plays a vital role in deriving error estimates of the control, state and co-state variables. Its proof is similar to that in [38,Theorem 4.1].

24)
where C > 0 is independent of h.
Then, using the approximation property of γ together with Lemmas 4 and 5 implies Adding (3.27) and (3.28) and using that ( y h − y h (u), γ ( y h − y h (u))) 0,Ω ≥ 0, we arrive at where we have used relations (2.2), (2.7), (2.11) and (2.13). An application of Lemmas 4 and 5 in the above inequality leads to the following bound Remark 1 (Right-hand side regularity) According to the contributions [8,19,28,30] (see also the references therein), for linear finite volume element methods applied to second order elliptic problems, the optimal error estimates (establishing second order accuracy in the L 2 − norm) can be achieved under the assumption that the source term is in H 1 (Ω). However, assuming that the right-hand side is in H 1 (Ω) does not imply that the exact solution is in H 3 (Ω), as discussed in e.g. [19]. Some counterexamples are actually given in [28,30] to confirm that optimal L 2 − error estimates cannot be derived if one only assumes that the forcing term is in L 2 (Ω). Proceeding analogously to the analysis of standard finite volume methods, optimal error estimates in the L 2 − norm have been derived by taking the source term in H 1 (Ω) (see for instance [19,24] and their references, for the specific case of DFV methods applied to elliptic and Stokes problems). Following the analysis of [31], one can derive the error estimates given in Lemmas 5 and 6 under the less restrictive assumption that f and y d are in H 1 (T ), that is, locally-H 1 .

Error Estimates Under Variational Discretisation
Theorem 1 Let ( y h , w h ) be DFV approximations of ( y, w) and let u h denote a variational discretisation of u. Then there exists a positive constant C independent of h, but depending on λ, such that the following estimates hold: Proof We recall the continuous variational inequality and the discrete variational inequality under variational discretisation and rearranging terms, we get

L 2 -Error Estimates for Fully Discretised Controls
A discrete admissible controlũ h = (ũ h, j ) d j=1 ∈ U h,ad is defined component-wise and locally asũ whereĨ h u j is the Lagrange interpolant of u j . To avoid ambiguity, we choose h sufficiently small so that min x∈T u j (x) = a j and max x∈T u j (x) = b j do not occur simultaneously within the same element T ∈ T h . Next, we proceed to group the elements in the primal mesh into three categories: h,n = ∅ for m = n according to the value of u j (x) on T . These sets are defined as On the other hand, the following assumption will be instrumental in the subsequent analysis. There exists a positive constant C independent of h such that (3.38) A similar assumption has been employed in [41][42][43]48]. We will first focus on error bounds for the control field under piecewise linear discretisation. Before proceeding we state an auxiliary result, whose proof can be found in [41]. (3.38) and that w ∈ W 1,∞ (Ω). Then, there exists C > 0 independent of h such that

Lemma 7 Assume
The main result in this section is stated as follows. Proof Testing the continuous and discrete variational inequalities against u h ∈ U h,ad ⊂ U ad andũ h ∈ U h,ad , respectively, and adding them, leads to Addition and subtraction ofũ h in the first term above yields and after rearranging terms we obtain In view of estimating the term u −ũ h 0,Ω , we proceed to rewrite it as whereas for T 2 , we employ the projection property (1.11) together with (3.38) Inserting the bounds of T 1 and T 2 in (3.40) we arrive at Finally, applying Cauchy-Schwarz and Young's inequalities, the estimates (3.24), (3.41), and Lemmas 5 and 7 into (3.39), we readily obtain the required result.
We now turn to the L 2 −error analysis for the control field under element-wise constant discretisation. The main idea follows from [14], using an L 2 −projection Π 0 : L 2 (Ω) −→ U h,0 that has the following property: there exists a positive constant C independent of h such that u − Π 0 u 0,Ω ≤ Ch u 1,Ω . Proof Since Π 0 U ad ⊂ U h,ad , the continuous and discrete optimality conditions readily imply that Adding and subtracting u, and rearranging terms, we then obtain and since Π 0 is an orthogonal projection and u h ∈ U h,ad , then the term λ(u h , Π 0 u − u) 0,Ω vanishes to give For the first term, we use (3.24) to get whereas a bound for I 2 follows from the orthogonality of Π 0 : It is left to show that w h is uniformly bounded, which can be readily derived using the coercivity of A h (·, ·) and c h (·, ·) and the uniform boundedness of U h,ad : Substituting the bounds for I 1 and I 2 in (3.43), and using (3.42) the desired result follows.

L 2 -Error Estimates for Velocity Under Full Discretisation of Control
The main result in this section is given as follows (see similar ideas, based on duality arguments also applied in [43,50]). ( y, w) be the state and co-state velocities, solutions of (1.1)-(1.4), and let ( y h , w h ) be their DFV approximations under piecewise linear (or piecewise constant) discretisation of control. Then

Theorem 4 Let
Proof We start by splitting the total error and applying triangle inequality as: where Π h represents the L 2 −projection operator onto the discrete control space U h . Next, let (w h ,r h ) ∈ V h × Q h be the unique solution of the auxiliary discrete dual Brinkman problem We then choosez h =w h andψ h =r h in (3.45) and (3.46), respectively, next we add the result, and we use the coercivity properties (2.8) and (2.12), to derive that (Π h u), respectively, and adding the result, we obtain (3.48) In addition, employing the discrete state equation for y h (u) and y h (Π h u), we obtain We then proceed to subtract (3.49) from (3.48) and to rearrange terms, to arrive at Using the definition of the norm |||·||| 0,h and its equivalence with the norm · 0,Ω we find that By virtue of the properties of Π h applied in the above inequality, we can assert that Approximation properties of γ and the L 2 −projection readily yield appropriate bounds for S 1 and S 2 , respectively: Then, a direct application of (3.47) yields We next use relations (2.11), (2.13) and (3.47) to obtain Consequently, substituting the estimates for the terms S 1 , S 2 , S 3 and S 4 back into (3.50), one straightforwardly arrives at The third term in (3.44) is bounded using (2.7) and proceeding as in the proof of Lemma 4: Using the discrete variational inequality along with the projection property of Π h and (3.37), we have the following relation Plugging the bounds for J 1 , J 2 , J 3 and J 4 in (3.53), putting (3.51) and (3.52) into (3.44), and using interpolation estimate u − Π h u 0,Ω ≤ Ch u 1,Ω along with Lemma 5; we can assert that y − y h 0,Ω ≤ Ch 2 y 2,Ω + p 1,Ω + u 1,Ω + f 1,Ω . (3.54) Finally, splitting the co-state velocity error as w − w h = w − w h ( y) + w h (y) − w h , using triangle inequality and Lemmas 4,5, and relation (3.54), we get the second desired estimate

Error Bounds in the Energy Norm
Theorem 5 Let ( y, w, p, r ) be the state and co-state velocities, and pressures, solutions of (1.1)-(1.4), and let ( y h , w h , p h , r h ) be their DFV approximations. Then Proof Using (3.5) and (3.6), applying triangle inequality and Lemma 4, we obtain The proof follows after combining Lemma 5 with the bounds for u − u h 0,Ω and y − y h 0,Ω .

Numerical Experiments
In this section we present a set of numerical examples to illustrate the theoretical results previously described. For sake of completeness, before jumping into the tests we provide some details about the implementation and algorithms for the efficient numerical realisation of the DFV method applied to the optimal control of Brinkman equations. Implementation Aspects We will use the well-known active set strategy (proposed in [6]) involving primal and dual variables (see also [25,44] for its application in Stokes flow). The principle is to approximate the constrained optimal control problem by a sequence of unconstrained problems, using active sets as summarised in Algorithm 1. By u n h , w n h we will denote the optimal control and adjoint velocity, solutions to the discrete problem (2.16)- (2.20) at the current iteration. Also, the control constraints are a = (a 1 , ..a d ) T and b = (b 1 , ..b d ) T . From u n h and w n h , construct the new finite sets A a n+1 , A b n+1 and I n+1 using (4.2) 5: if n ≥ 1, and A a n+1 = A a n , and where χ T * i is the characteristic function assuming the value 1 on T * i ∈ T * h and zero elsewhere. We next proceed to define the discrete active and inactive sets, based on the degrees of freedom of U h , as follows where, in general, s n,k j,h stands for the discrete value associated with the degree of freedom at position k, related to the spatial component j of the vector field s, at the step n of Algorithm 1. By the definition of the optimal control problem, we have that and if we further introduce the following characteristic sets (4.3)  Finally, we define the matrix blocks along with the vectors so that after testing (4.3) against {ψ i } M i=1 we end up with the following matrix form of the discrete optimal control problem (2.16)- (2.20): where Y, P, W, R and U are the coefficients in the expansion of y n+1 h , p n+1 h , w n+1 h , r n+1 h and u n+1 h , respectively, and the hats indicate quantities associated with the previous iteration.
Example 1 We start by assessing the experimental convergence of the proposed scheme applied to the optimal control problem (1.1)-(1.4) defined on the unit square Ω = (0, 1) 2 . Viscosity, permeability and the weight for the control cost assume the following constant (see e.g. [48]) which satisfy the homogeneous Dirichlet boundary conditions under which the analysis was performed. Source term and desired velocity field of the problem are constructed according to these exact solutions, that is, respectively A family of nested primal and dual triangulations of Ω is generated, on which we compute errors in the L 2 −and mesh-dependent norm |||·||| 1,h for the state and co-state velocity, in the L 2 −norm for pressures, and in the L 2 −norm for the control approximation. Table 1 displays the error history for this first test, where we observe optimal convergence rates for velocity and pressure (only those of the state equation are shown) in their natural norms, along with an O(h) convergence for the control when approximated by piecewise constant elements, which improves to roughly a O(h 3/2 ) rate under piecewise linear approximations. We can also confirm that a maximum of three iterations are needed to reach the stopping criterion that the active sets are equal to those in the previous optimisation step. This indicates a mesh independence of the method in the sense that the number of iterations needed to achieve the stopping criterion is independent of the resolution. In addition we portray in Fig. 2 the obtained approximate solutions at the finest resolution level, where we highlight the active sets with a contour plot on top of the control and state velocities. In all examples herein we employ a BiCGSTAB method with AMG preconditioning to solve the linear systems involved at each step of Algorithm 1. Moreover, the zero-mean pressure condition is applied for both pressure and adjoint pressure using a real Lagrange multiplier approach (which accounts for adding one column and one row to the relevant matrix system).
We also present a basic comparison with other classical methods in terms of accuracy. For instance, we have performed the same test as above but employing model coefficients with jumps, in order to highlight the need for discontinuous approximations. Both fluid viscosity and medium permeability have now a discontinuity of five orders of magnitude at x 1 = 0.5. The tested methods are: a conforming stable P 2 −P 0 and MINI-element pairs for velocity and pressure approximation, a classical interior penalty DG method using the same stabilisation parameters as in (2.4)-(2.5), and the proposed DFV formulation. In all cases we consider a piecewise linear approximation of the control variable. The results are collected in Fig. 3, where convergence histories (errors for velocity and pressure vs. the number of degrees of freedom DoF= 2(N + L) + M) associated to the studied discretisations are shown. For all fields, the DFV approximation exhibits a slightly better accuracy than its pure-DG counterpart. This may be explained by the smaller elements used in the dual mesh (but being associated with the same number of DoF). On the other hand, for coarse meshes the conforming approximation P 2 − P 0 outperforms all other methods, but for finer meshes the discontinuous coefficients of the problem imply a badly conditioned system matrix requiring more iterations of the linear solver and eventually the conforming methods lose their optimal convergence. For a fixed number of DoF, the proposed DFV scheme produces smaller errors for the pressure approximation than the other methods. We stress that some recent theoretical comparison results are available for forward Stokes problems (see e.g. [13]), but only in the case of smooth solutions and constant coefficients. If the comparisons are carried out for the case of smooth solutions, then the error estimates in e.g. Theorem 5 are indeed of the same order as their finite element counterpart. However the constants in the estimates are not necessarily the same. As mentioned above, since the dual elements are in principle smaller than the primal ones, the approximate solutions and the corresponding errors generated with discontinuous finite volume schemes are still slightly more localised, implying that the errors themselves are smaller than those produced with methods based on the primal mesh. Complexity, implementation, and CPU times for assembly  and solution of the linear systems are, on the other hand, comparable to those associated with classical finite elements.
Example 2 Our second test focuses on the optimal control problem applied to the well-known lid driven cavity problem. The objective function still corresponds to (1.1), but no analytic exact solution is available. Again the domain consists of the unit square, and the data of the problem are given by a traction boundary condition on the top of the lid, the applied body force, and an observed velocity field y = (1, 0) T on the top and zero elsewhere, and f = y d = 0 in Ω. The adjoint problem is subject to homogeneous Dirichlet data. The viscosity is set to ν = 0.1, the control weight is now λ = 0.2, the admissible control space is characterised by a 1 = a 2 = −0.15, b 1 = b 2 = 0.15, and the permeability exhibits a discontinuity on the line x 2 = 0.4: K(x) = κ ν I, with κ(x) = {10,000 if x 2 ≥ 0.4; 10 elsewhere in Ω} (see also [2,3] for the simulation of Brinkman flows with sharp interfaces). The domain is discretised into 20,000 primal triangular elements, and Fig. 4 portrays all fields obtained with our DFV scheme, where the stabilisation parameter is α d = 10. From Fig. 4 we observe that the controlled velocity approaches the desired velocity, that is, it goes to zero and the movement of the fluid concentrates in the upper section of the cavity. In addition, we study the influence of the Tikhonov regularisation in the iteration count of the active set algorithm applied to a coarse solve of this test. As in [25], we immediately observe that a larger number of iterations are required for smaller values of λ (see Table 2).
Example 3 Next we turn to the numerical solution of a three-dimensional optimal control problem (see also [34]). The domain is a cylinder with height 4 and radius 1, aligned with the x 2 axis. The permeability field is anisotropic K = diag(0.1, 10 −5 χ B + 0.1χ B c , 0.1), where B is a ball of radius 1/4 located at the centre of the domain. A Poiseuille inflow profile is imposed as state velocity at x 2 = 0: y = (0, 10(1−x 1 2 −(x 3 −1/2) 2 ), 0) T , a zero-pressure is considered on x 2 = 4, whereas homogeneous Dirichlet data are enforced on the remainder of ∂Ω. The viscosity is ν = 0.005, the Tikhonov regularisation is λ = 1/2, the desired velocity is zero y d = 0, the bounds for the control are a j = a = −0.1 and b j = b = 0.2, and a smooth body force is set as in [4]: f = K −1 (exp(−x 2 x 3 ) + x 1 exp (−x 2 2 ), cos(π x 1 ) cos(π x 3 ) − x 2 exp (−x 2 2 ), −x 1 x 2 x 3 − x 3 exp(−x 2 3 )) T . The primal meshes has 78,631 internal tetrahedral elements and 13,593 vertices. We observe that five iterations are required to reach the stopping criterion (4.1). Snapshots of the resulting approximate fields are collected in Fig. 5.