SUPG-Stabilized Virtual Element Method for Optimal Control Problem Governed by a Convection Dominated Diffusion Equation

In this paper, the streamline upwind/Petrov Galerkin (SUPG) stabilized virtual element method (VEM) for optimal control problem governed by a convection dominated diffusion equation is investigated. The virtual element discrete scheme is constructed based on the first-optimize-then-discretize strategy and SUPG stabilized virtual element approximation of the state equation and adjoint state equation. An a priori error estimate is derived for both the state, adjoint state, and the control. Numerical experiments are carried out to illustrate the theoretical findings.


Introduction
In this paper, we aim to discuss a priori error analysis of SUPG stabilized virtual element method (VEM) for the optimal control problem governed by the convection dominated diffusion equation. We consider the following optimal control problem: subject to −∇ · (ε∇y) + β(x) · ∇y + δy = f + u in Ω, where J(y, u) is the objective functional, y d ∈ L 2 (Ω) is the desired state, and γ > 0 is the regularization parameter. 0 < ε 1 represents constant diffusion coefficient and β ∈ [W 1,∞ (Ω)] 2 with ∇ · β = 0 is the transport advective field. δ > 0 is a constant. f ∈ L 2 (Ω) is the volume source term. Ω ⊂ R 2 is a bounded domain with Γ = ∂Ω. The control constraint set is given by U ad = {u ∈ L 2 (Ω) : u a ≤ u(x) ≤ u b a.e. in Ω with u a , u b ∈ R and u a ≤ u b }.
Optimal control problems governed by the convection dominated diffusion equation have many applications in real life, such as the air pollution problem ( [1]) and waste water treatment problem ( [2]). It is well known that a characteristic of convection-dominated equation is that the solutions may have the sharp boundary and interior layers when the coefficient of the convection field is relatively large. Since numerical methods without any treatment do not work well in this case, various robust schemes such as SUPG formulation, residual-free bubbles methods, and discontinuous Galerkin methods for convective dominance equations have been developed. For the numerical approximation of the convection diffusion optimal control problem, we refer to [3][4][5][6][7][8][9][10].
VEM can be regarded as a generalization of the finite element method (FEM) to general polygonal and polyhedral meshes, and it is originally introduced in [11] as a C 0 -conforming method for solving the two-dimensional Poisson equation. Thus far, VEM has been used in a variety of problems, such as elliptic problem ( [12][13][14]), parabolic problem ( [15]), and Stokes problem ( [16,17]). The method performs well in geometrically complex domains [18] and with badly shaped polygonal elements [19]. In fact, the underlying virtual element space can be seen as the finite element space plus some suitable non-polynomial functions, which are the solutions of PDE problems inside each element. Compared with SUPG-FEM, SUPG-VEM provides great flexibility for us to use arbitrary polygonal meshes (even nonconvex). SUPG-VEM also has great advantages for adaptive refinement. For instance, locally adapted meshes do not require any local post processing because polygonal meshes are allowed, and any limitations caused by maximum angle conditions or mesh distortion are eliminated. In addition, we do not need to add additional degrees of freedom for hanging nodes during adaptive refinement, since we can just treat the hanging nodes as new nodes. However, there are not many studies on SUPG-VEM of convection dominant problems. Cangiani et al. first studied the non-consistent SUPG-VEM problem of convection-dominated diffusion in [20]. Subsequently, SUPG-stabilized conforming and non-conforming VEMs are presented in [21,22]. However, the stability and convergence analysis in [21,22] is not uniform in the diffusion/convection parameters, and a small enough mesh size is needed for analysis. Recently, Beirão da Veiga et al. discussed a robustness analysis of SUPG-stabilized virtual elements for diffusion-convection problems in [23]. By slightly modifying the SUPG format of [21], they propose a new way to discretize the convection term, which ultimately demonstrates the robustness of the parameters involved in the convergence estimation without requiring sufficiently small mesh sizes.
As we know the application of the virtual element method in the optimal control problem was not reported up to now. Compared with FEM, the computability of the discrete scheme is more important since the virtual element space contains non-polynominal functions. By projection operators and the first-optimize-then-discretize strategy, we construct a computable SUPG-stabilized VEM discrete scheme for the optimal control problem governed by the convection dominated diffusion equation, where the control is implicitly discretized. Moreover, inspired by [23], we use a novel discretization of the convection term that allows us to develop error estimates that are fully robust in the convection dominated cases. We derive an a priori error estimate for the optimal control problem by introducing some auxiliary problems and present a projected gradient algorithm to solve the discrete optimal control problem. Finally, we carry out some numerical examples to verify our theoretical analysis.
The paper is organized as follows. In the next section, we give some preliminary knowledge about virtual element space and the projection operators. In Section 3, the SUPG stabilizing virtual element discrete scheme is constructed. In Section 4, a priori error estimates are derived for the state, adjoint state, and control. In Section 5, we perform some numerical experiments to verify the theoretical results.
Throughout the paper, the symbol denotes a bound up to a generic positive constant, independent of the mesh size h, of the SUPG parameter τ E , of the diffusive coefficient ε and of the transport advective field β. Moreover, the analysis of the three-dimensional case could be developed with very similar arguments.

Preliminaries
In this section, we introduce the virtual element space and projection operators. Let T h be a family of decompositions of the domain Ω into star-shaped polygonals E. h E denotes the diameter of element E, i.e., the maximum distance between any two points on element E and h = sup E∈T h h E . We further assume that ∂E denotes the edges of E ∈T h . n E and h s denote the unit outward normal vector to ∂E and the length of edge s, respectively. The following assumption about the mesh is necessary for the theoretical analysis. Assumption 1 ([12]). We assume the existence of a constant ρ > 0 such that, for all h > 0 and for all E ∈ T h ,

•
Every element E of T h is star-shaped with respect to a ball of radius bigger or equal to ρh E ; • For every element E of T h and every edge s of E, h s ≥ ρh E .
For n = 0, 1, P n (E) denotes the space of polynomials of degree ≤ n on E (with P −∞ (E) = {0}) and the following polynomial projections [12] are given: with an obvious extension for vector functions Π 0 n : plus (to take care of the constant part of Π 1 n v): Then, following [12], the global virtual element space is defined as where the space P 1 (E)/P −∞ (E) denotes the polynomials in P 1 (E) that are L 2 (E)-orthogonal to P −∞ (E), i.e., the degree of the polynomials in this space are 1 and 0. Following [12], for each element E ∈ T h , the local virtual element space V E h contains the space P 1 (E). The vertices of a polygonal element E with N E edges are denoted by υ i for i = 1, ..., N E . Because the polygonal element E has no internal degrees of freedom, we can simply choose the υ i for i = 1, ..., N E to be the degrees of freedom of polygonal element E. Remark 1. Following [12], from the definition of the virtual element space and the two projection operators, we have Π 1 1 = Π 0 1 .

Lemma 1 ([23]
). Under the Assumption 1, for any E ∈ T h and for any smooth enough function ϕ defined on E, it holds

SUPG Stabilizing Virtual Element Approximation for Optimal Control Problem
For all y, v ∈ H 1 0 (Ω), we set By a direct computation, being ∇ · β = 0, it is easy to see that the bilinear form b(·, ·) is skew symmetric, i.e., Therefore, we can rewrite the bilinear form b(·, ·) as Then, the weak form of optimal control Equations (1) and (2) is characterized by Here, the convection term is rewritten in a skew symmetric form, which is a useful step for the VEM to ensure that the discrete framework preserves the properties of the symmetric and skew-symmetric parts of the bilinear form, see [14,23].
To derive the first order optimality system, we define the Lagrangian functional as follows: L(y, p, u) = J(y, u) + ( f + u, p) − A(y, p).
Taking the directional derivative for L(y, p, u) with respect to y, p, and u, we obtain the continuous first order optimality system of Equations (1) and (2) as follows: denote the pointwise projection onto the admissible set U ad . The optimal inequality is equivalent to u = P U ad (− 1 γ p). We can derive that the adjoint state equation has the strong form in Ω, The virtual element discrete scheme of state equation with SUPG stabilization can be defined as follows: find y h (u) ∈ V h , such that Here, Following [11], S E is any symmetric positive definite bilinear form to be chosen to verify where α * and α * are two positive constants independent of E and h E . There are many choices for S E , and, following [11], we take the simple choice where dof r (y h ) denotes the value of the rth local degree of freedom defining y h in V E h .

Remark 2.
Due to limited regularity of the control problem, we restrict k = 1 in the virtual element space. In this case, we can observe that Q E h (·, ·) = 0 by the definition of the L 2 projection operator.
Similar to the state equation, we also use SUPG-stabilized VEM to discretize the adjoint state equation. Then, we can define the discrete first order optimality system as follows: where

A Priori Error Estimates
In this section, we first define the VEM SUPG norm and introduce the auxiliary problems. Then, under certain data assumption, the error estimate of the auxiliary problem and the optimal control problem in the VEM SUPG norm are given. Finally, we derive the error estimate between Equations (4) and (6).
We first define the local VEM SUPG norm Notice that the norm · supg,E here is slightly different from the standard SUPG norm · SUPG,E introduced in standard SUPG theory, i.e., However, for all w h ∈ H 1 (E), using the fact that This implies that the standard SUPG norm can be controlled by the VEM SUPG norm.
To derive an a priori error estimate, we need to introduce the following auxiliary problems: and Assumption 2 (Data assumption). The solutions y, p, u of the optimal control problem and the f , y d satisfy: Note that (y h (u), p h (y)) are SUPG VEM approximation of (y, p). The following results are not restrictive to β E > 0 since β E = 0 implies β| E = 0 and thus the corresponding terms vanish.

Lemma 3 ([23]
). Let (y, p) and (y h (u), p h (y)) be the solutions of (4) and (9), respectively. Then, the following estimates hold under the Assumptions 1, 2 and, in the case of a convection dominated Next, we derive a priori error estimates for SUPG VEM approximation of the optimal control problem. Theorem 1. Let (y, p, u) and (y h , p h , u h ) be the solutions of (4) and (6), respectively. Then, we have that Note that By the auxiliary Equation (8), we deduce Then, we have This shows Note that Then, by Lemma 1, we have

Theorem 2.
Under the Assumptions 1 and 2 and, in case of a convection dominated regime, let (y, p, u) and (y h , p h , u h ) be the solutions of Equations (4) and (6), respectively. Then, we have following estimate: Proof. We decompose the errors y − y h and p − p h into By the Lemma 3, we have Moreover, by the governing equation of y h , y h (u) in (6) and (9), respectively, we have Recalling the coercivity of A supg,h in (7), Π 0 1 = Π 1 1 and the fact that τ E ≤ C τ δ we can get This implies y h (u) − y h supg u − u h 0,Ω .
Combining the above inequalities gives By Lemma 3, we also have Similar to the estimate (10), by the definition of p h (y), p h (u) in (9), we also derive This implies From the coercivity of B and definition of p(y h ), p(u h ) in (8), we can deduce Thus, we can get p(y h ) − p(u h ) 0,Ω ≤ y h − y(u h ) 0,Ω . Using the results of [23], we obtain By the triangle inequality and the relationship between · SUPG and · supg , we arrive at From Theorem 1 and Equation (11), we have the following estimate: Inserting the above estimate into the estimates of state and adjoint state yields the final result.

Numerical Results
In this section, we give the mesh types of our numerical experiments and introduce a projected gradient algorithm based on the SUPG stabilized discrete first order optimality system (6) to verify our a priori analysis. Since the VEM solutions y h and p h are not explicitly known inside the elements, we use e y,0 , e y,1 , c y to represent the L 2 norm, H 1 norm, and standard SUPG norm between y and Π 1 1 y h . Similarly, e p,0 , e p,1 , c p denote the L 2 norm, H 1 norm and standard SUPG norm between p and Π 1 1 p h , and e u,0 denotes the L 2 norm between u and u h .
For the mesh types, we consider distorted square, hexagonal, and Lloyd mesh ( [24]), which are shown as Figure 1a The projected gradient algorithm is given below (Algorithm 1): Algorithm 1 Projected gradient algorithm. Require: Regularization parameter γ and tolerance error η. Ensure: 1.
Choose the initial value u h . Set error = 1.
Solving the state equation in the discrete first order optimality system (6) to get state variable y h ; 4.
Solving the adjoint state equation in the discrete first order optimality system (6) to obtain adjoint state variable p h ; 5.
Control variable u hnew are obtained by using the projection P U ad ;
end while For the detailed calculation process of the projection and the influence of projection on the convergence rate, we can refer to the literature [14]. (1) and (2)

Example 1. Consider the optimal control Equations
The exact solutions are chosen to be f and y d can be determined from the exact solutions y, p, u.
We choose C τ = 0.8 and let δ under a fixed mesh. In order to ensure the dominance of the convection term, i.e., the choice of parameter , numerical comparisons of r 1 , r 2 , and r 3 are shown in Tables 1 and 2 for each type of grid subdivision with ε = 10 −10 . We can see that the convection term is always dominant. We also give the convergence rates of the above norms for distorted square, hexagonal, and Lloyd mesh with the SUPG term. We can observe that the convergence orders of L 2 errors, H 1 errors, and standard convective norm errors approximate 2, 1, and 1.5, respectively, which are the optimal convergence rates. The experimental data verify our theoretical analysis. In Figure 2a-c, we plot the profiles of the exact solutions of state, adjoint state, and control with ε = 10 −10 , respectively. Figure 2 shows an intuitive comparison between the unstabilized numerically computed state, adjoint state, control and the numerically computed solutions obtained using the SUPG stabilization on Lloyd mesh with ε = 10 −10 . By comparison, we can find that the numerical solutions Figure 2d-f show a very good agreement with the exact solutions and show a good stability as ε → 0 when the SUPG term exists. The quality of the discrete solutions deteriorates obviously when there is no SUPG term (Figure 2g-i).  . The exact solution of the optimal control problem is as follows: These functions are inserted into the equations and the corresponding source terms f , y 0 , and u 0 are computed.
We also choose C τ = 0.8 and let δ under a fixed mesh. In Figure 3, we show the convergence graphs for three meshes with an SUPG term under the above norms and the intuitive comparison of r 1 , r 2 and r 3 , respectively. We can observe that the convergence orders of L 2 errors, H 1 errors and standard SUPG norm errors are approximately parallel to the lines with slopes 2, 1, and 1.5 in the convection dominated regime. We can observe that the convergence rates are in agreement with the theoretical prediction.  The contour-lines and profiles of the stabilized numerically computed state, control, and the profiles of the numerically computed state and control without SUPG term on Lloyd mesh are shown in Figure 4, respectively. We can observe that the quality of the numerical solutions are good when the SUPG term is present; otherwise, the numerical solutions are obviously destroyed, which implies that the SUPG term has a good effect.
The right-hand term f and the desired state y d can be calculated by the exact solutions and governing equations.
In this example, the state y and adjoint state p have sharp interior layers along the circles (x 1 − 1) 2 + (x 2 − 1) 2 = 0.7 2 and x 2 1 + x 2 2 = 0.5 2 , respectively. Figure 5 shows the contour-lines and profiles of stabilized numerical solutions on the distorted square mesh and the quality of the numerical solutions are not destroyed by the numerical oscillation. We can see that our method represents and processes the interior layers well. The corresponding convergence rates are given in Figure 6, which are also consistent with our theoretical analysis.

Conclusions
In this paper, we attempt to apply SUPG-stabilized VEM to approximate an optimal control problem governed by a convection dominated diffusion equation with pointwise control constraint. A priori error estimates are derived. The theoretical findings are verified by numerical examples.
Since the VEM has great flexibility in the mesh partition, in our future work, we are going to investigate VEM approximation of an optimal control problem governed by fractional advection-diffusion-reaction equations ( [25][26][27]).