Next Article in Journal
Water-Spray-Cooled Quasi-Isothermal Compression Method: Water-Spray Flow Improvement
Next Article in Special Issue
Existence of a Unique Solution to a Fractional Partial Differential Equation and Its Continuous Dependence on Parameters
Previous Article in Journal
Minimax Rates of p-Losses for High-Dimensional Linear Errors-in-Variables Models over q-Balls
Previous Article in Special Issue
Stability of Non-Linear Dirichlet Problems with ϕ-Laplacian
Article

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

School of Mathematics and Statistics, Shandong Normal University, Jinan 250000, China
*
Author to whom correspondence should be addressed.
Academic Editors: Ravi P. Agarwal and Maria Alessandra Ragusa
Received: 5 April 2021 / Revised: 24 May 2021 / Accepted: 30 May 2021 / Published: 5 June 2021
(This article belongs to the Special Issue Nonlinear Dynamics and Analysis)

Abstract

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.
Keywords: optimal control; convection dominated diffusion equation; SUPG stabilized VEM; a priori error estimate optimal control; convection dominated diffusion equation; SUPG stabilized VEM; a priori error estimate

1. 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:
min u U a d J ( y , u ) : = 1 2 y y d 2 + γ 2 u 2
subject to
· ( ε y ) + β ( x ) · y + δ y = f + u in Ω , y = 0 on Γ ,
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 a d = { 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 non-convex). 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.

2. 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:
  • the L 2 p r o j e c t i o n Π n 0 : L 2 ( E ) P n ( E ) , defined by
    ( q , v ) 0 , E = ( q , Π n 0 v ) 0 , E f o r a l l v L 2 ( E ) a n d q P n ( E ) ,
    with an obvious extension for vector functions Π n 0 : [ L 2 ( E ) ] 2 [ P n ( E ) ] 2 ;
  • the H 1 p r o j e c t i o n Π n 1 : H 1 ( E ) P n ( E ) , defined by
    ( q , v ) 0 , E = ( q , Π n 1 v ) 0 , E f o r a l l v H 1 ( E ) a n d q P n ( E )
    plus (to take care of the constant part of Π n 1 v ):
    E ( v Π n 1 v ) d s = 0 .
Then, following [12], the global virtual element space is defined as
V h : = { v h H 0 1 ( Ω ) : v h | E V h E , E T h } ,
where
V h E : = { v h H 1 ( E ) C 0 ( E ) : Δ v h E P 1 ( E ) , v h | E P 1 ( s ) s E , ( v h , p ) 0 , E = ( Π 1 1 v h , p ) 0 , E p P 1 ( E ) / P ( E ) } ,
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 h E 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 = Π 1 0 .
Lemma 1
([23]). Under the Assumption 1, for any E T h and for any smooth enough function φ defined on E, it holds
φ Π n 0 φ W p m ( E ) h E s m | φ | W p s ( E ) s , m N , m s n + 1 , p = 1 , . . . , , φ Π n 1 φ m , E h E s m | φ | s , E s , m N , m s n + 1 , s 1 , φ Π n 0 φ m , E h E s 1 m | φ | s , E s , m N , m s n + 2 , s 1 .

3. SUPG Stabilizing Virtual Element Approximation for Optimal Control Problem

For all y , v H 0 1 ( Ω ) , we set
a ( y , v ) : = Ω ε y · v d Ω , b ( y , v ) : = Ω β · y v d Ω , c ( y , v ) : = Ω δ y v d Ω .
By a direct computation, being · β = 0 , it is easy to see that the bilinear form b ( · , · ) is skew symmetric, i.e.,
b ( y , v ) = b ( v , y ) f o r a l l y , v H 0 1 ( Ω ) .
Therefore, we can rewrite the bilinear form b ( · , · ) as
b ( y , v ) = 1 2 b ( y , v ) b ( v , y ) .
Then, the weak form of optimal control Equations (1) and (2) is characterized by
min u U a d J ( y , u ) : = 1 2 y y d 2 + γ 2 u 2
subject to
A ( y , v ) = ( f + u , v ) v H 0 1 ( Ω ) ,
where
A ( y , v ) : = a ( y , v ) + b ( y , v ) + c ( y , v ) .
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:
A ( y , w ) = ( f + u , w ) , w H 0 1 ( Ω ) , B ( p , w ) = ( y y d , w ) , w H 0 1 ( Ω ) , ( γ u + p , v u ) 0 , v U a d ,
where B ( p , w ) : = a ( p , w ) b ( p , w ) + c ( p , w ) .
Let
P U a d ( u ) = max { u a , min { u , u b } }
denote the pointwise projection onto the admissible set U a d . The optimal inequality is equivalent to u = P U a d ( 1 γ p ) . We can derive that the adjoint state equation has the strong form
· ( ε p ) β ( x ) · p + δ p = y y d in Ω , p = 0 on Γ .
The virtual element discrete scheme of state equation with SUPG stabilization can be defined as follows: find y h ( u ) V h , such that
A supg , h ( y h ( u ) , w h ) = E T h ( Π 1 0 ( f + u ) , w h + τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h .
Here,
A supg , h ( v h , w h ) : = E T h A supg , h E ( v h , w h ) : = E T h ( a h E ( v h , w h ) + b h E ( v h , w h ) + c h E ( v h , w h ) + Q h E ( v h , w h ) + B h E ( v h , w h ) + R h E ( v h , w h ) ) ,
where
a h E ( v h , w h ) : = E ε Π 0 0 v h · Π 0 0 w h d E + ε S E ( ( I Π 1 1 ) v h , ( I Π 1 1 ) w h ) , b h E ( v h , w h ) : = 1 2 [ E β · Π 1 0 v h Π 1 0 w h d E E β · Π 1 0 w h Π 1 0 v h d E + E ( β · n E ) ( I Π 1 0 ) v h Π 1 0 w h d s E ( β · n E ) ( I Π 1 0 ) w h Π 1 0 v h d s ] , c h E ( v h , w h ) : = E δ Π 1 0 v h Π 1 0 w h d E , Q h E ( v h , w h ) : = τ E E ε · Π 0 0 v h ( β · Π 0 0 w h ) d E , B h E ( v h , w h ) : = τ E E β · Π 0 0 v h β · Π 0 0 w h d E + τ E β E 2 S E ( ( I Π 1 1 ) v h , ( I Π 1 1 ) w h ) , R h E ( v h , w h ) : = τ E E δ Π 1 0 v h ( β · Π 0 0 w h ) d E ,
τ E > 0 is the SUPG parameter and β E = β [ L ( E ) ] 2 . Following [11], S E is any symmetric positive definite bilinear form to be chosen to verify
α | w h | 1 , E 2 S E ( w h , w h ) α | w h | 1 , E 2 f o r a l l w h K e r ( Π 1 1 ) ,
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
S E ( y h , w h ) = r = 1 N E dof r ( y h ) dof r ( w h ) ,
where dof r ( y h ) denotes the value of the rth local degree of freedom defining y h in V h E .
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 h E ( · , · ) = 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: find ( u h , y h , p h ) U a d × V h × V h , such that
A supg , h ( y h , w h ) = E T h ( Π 1 0 ( f + u h ) , w h + τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h , B supg , h ( p h , w h ) = E T h ( Π 1 0 ( y h y d ) , w h τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h , E T h ( γ u h + Π 1 0 p h , v h u h ) 0 , E 0 , v h U a d ,
where
B supg , h ( p h , w h ) : = E T h B supg , h E ( p h , w h ) : = E T h ( a h E ( p h , w h ) b h E ( p h , w h ) + c h E ( p h , w h ) Q h E ( p h , w h ) + B h E ( p h , w h ) R h E ( p h , w h ) ) .

4. 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
w h supg , E 2 : = ε w h 0 , E 2 + τ E β · Π 0 0 w h 0 , E 2 + τ E β E 2 ( I Π 1 1 ) w h 0 , E 2 + δ ( I Π 1 1 ) w h 0 , E 2 + δ Π 1 1 w h 0 , E 2
for all w h H 1 ( E ) and the global VEM SUPG norm w h supg 2 : = E T h w h supg , E 2 for all w h H 1 ( Ω ) . Notice that the norm · supg , E here is slightly different from the standard SUPG norm · SUPG , E introduced in standard SUPG theory, i.e.,
w h SUPG , E 2 : = ε w h 0 , E 2 + τ E β · w h 0 , E 2 + δ w h 0 , E 2 .
However, for all w h H 1 ( E ) , using the fact that
p h Π 0 0 p h 0 , E = ( I Π 0 0 ) ( p h Π 1 1 p h ) 0 , E p h Π 1 1 p h 0 , E
we arrive at
β · w h 0 , E 2 = β · ( w h Π 0 0 w h + Π 0 0 w h ) 0 , E 2 2 β · Π 0 0 w h 0 , E 2 + 2 β E 2 ( I Π 0 0 ) w h 0 , E 2 2 β · Π 0 0 w h 0 , E 2 + 2 β E 2 ( I Π 1 1 ) w h 0 , E 2
and
δ w h 0 , E 2 δ ( w h Π 1 1 w h + Π 1 1 w h ) 0 , E 2 2 δ Π 1 1 w h 0 , E 2 + 2 δ ( I Π 1 1 ) w h 0 , E 2 .
This implies that the standard SUPG norm can be controlled by the VEM SUPG norm.
Lemma 2
([23]). Under the Assumption 1, if there exists a constant C τ ( 0 , 2 ) such that the parameters τ E satisfy τ E C τ δ , E T h , the bilinear form A supg , h E ( · , · ) and B supg , h E ( · , · ) satisfy for all w h V h ( E ) the coercivity inequality
w h supg , E 2 A supg , h E ( w h , w h ) , w h supg , E 2 B supg , h E ( w h , w h ) .
To derive an a priori error estimate, we need to introduce the following auxiliary problems:
A ( y ( u h ) , w ) = ( f + u h , w ) , w H 0 1 ( Ω ) , B ( p ( y h ) , w ) = ( y h y d , w ) , w H 0 1 ( Ω ) , B ( p ( u h ) , w ) = ( y ( u h ) y d , w ) , w H 0 1 ( Ω )
and
A supg , h ( y h ( u ) , w h ) = E T h ( Π 1 0 ( f + u ) , w h + τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h , B supg , h ( p h ( u ) , w h ) = E T h ( Π 1 0 ( y h ( u ) y d ) , w h τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h , B supg , h ( p h ( y ) , w h ) = E T h ( Π 1 0 ( y y d ) , w h τ E ( β · Π 0 0 w h ) ) 0 , E , w h V h .
Assumption 2
(Data assumption). The solutions y , p , u of the optimal control problem and the f , y d satisfy:
f , u , y d H 1 ( T h ) , y , p H 2 ( T h ) .
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 regime, i.e., τ E = min { h E β E , h E 2 ε , C τ δ } = h E β E
y y h ( u ) supg 2 + p p h ( y ) supg 2 E T h ( h E 3 ( β E + β E 1 ) + ε h E 2 + β E 1 ε 2 h E + h E 4 + β E 1 h E 5 ) E T h ( h E 3 ( β E + β E 1 + h E + β E 1 h E 2 ) ) = O ( h 3 ) .
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
u u h 0 , Ω h 2 + p h p ( u h ) 0 , Ω .
Proof. 
Define
J ^ ( u ) ( v u ) = Ω ( γ u + p ( u ) ) ( v u ) d x .
Note that
J ^ ( u ) ( u u h ) J ^ ( u h ) ( u u h ) = γ Ω ( u u h ) 2 d x + Ω ( p ( u ) p ( u h ) ) ( u u h ) d x .
By the auxiliary Equation (8), we deduce
Ω ( p ( u ) p ( u h ) ) ( u u h ) d x = A ( y ( u ) y ( u h ) , p ( u ) p ( u h ) ) = B ( p ( u ) p ( u h ) , y ( u ) y ( u h ) ) = ( y ( u ) y ( u h ) , y ( u ) y ( u h ) ) 0 , Ω 0 .
Then, we have
γ u u h 0 , Ω 2 J ^ ( u ) ( u u h ) J ^ ( u h ) ( u u h ) = Ω ( γ u + p γ u h p ( u h ) ) ( u u h ) d x = ( γ u + p , u u h ) 0 , Ω + E T h ( γ u h + Π 1 0 p h , u h u ) 0 , E + E T h ( Π 1 0 p h p ( u h ) , u u h ) 0 , E 0 + 0 + E T h ( Π 1 0 p h p ( u h ) , u u h ) 0 , E .
This shows
u u h 0 , Ω E T h Π 1 0 p h p ( u h ) 0 , E 2 1 2 .
Note that
Π 1 0 p h p ( u h ) 0 , E   Π 1 0 p h Π 1 0 p ( u h ) 0 , E + Π 1 0 p ( u h ) p ( u h ) 0 , E p h p ( u h ) 0 , E + Π 1 0 p ( u h ) p ( u h ) 0 , E .
Then, by Lemma 1, we have
u u h 0 , Ω h 2 + p h p ( u h ) 0 , Ω .
   □
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:
y y h supg + p p h supg + u u h 0 , Ω h 3 / 2 .
Proof. 
We decompose the errors y y h and p p h into
y y h = y y h ( u ) + y h ( u ) y h , p p h = p p h ( u ) + p h ( u ) p h .
By the Lemma 3, we have
y y h ( u ) supg h 3 / 2 .
Moreover, by the governing equation of y h , y h ( u ) in (6) and (9), respectively, we have
A supg , h ( y h ( u ) y h , w h ) = E T h ( Π 1 0 ( u u h ) , w h + τ E β · Π 0 0 w h ) 0 , E .
Let w h = y h ( u ) y h . Recalling the coercivity of A supg , h in (7), Π 1 0 = Π 1 1 and the fact that τ E C τ δ we can get
y h ( u ) y h supg 2 A supg , h ( y h ( u ) y h , y h ( u ) y h ) = E T h ( Π 1 0 ( u u h ) , y h ( u ) y h + τ E β · Π 0 0 ( y h ( u ) y h ) ) 0 , E u u h 0 , Ω ( E T h Π 1 0 ( y h ( u ) y h ) 0 , E 2 ) 1 / 2 + u u h 0 , Ω ( E T h τ E β · Π 0 0 ( y h ( u ) y h ) 0 , E 2 ) 1 / 2 u u h 0 , Ω y h ( u ) y h supg + u u h 0 , Ω ( E T h τ E y h ( u ) y h supg , E 2 ) 1 / 2 u u h 0 , Ω y h ( u ) y h supg + u u h 0 , Ω ( E T h C τ δ y h ( u ) y h supg , E 2 ) 1 / 2 u u h 0 , Ω y h ( u ) y h supg .
This implies
y h ( u ) y h supg u u h 0 , Ω .
Combining the above inequalities gives
y y h supg h 3 / 2 + u u h 0 , Ω .
By Lemma 3, we also have
p p h ( y ) supg h 3 / 2 .
Similar to the estimate (10), by the definition of p h ( y ) , p h ( u ) in (9), we also derive
p h ( y ) p h ( u ) supg y y h ( u ) supg h 3 / 2
and
p h ( u ) p h supg y h ( u ) y h supg u u h 0 , Ω .
This implies
p p h supg h 3 / 2 + u u h 0 , Ω .
From the coercivity of B and definition of p ( y h ) , p ( u h ) in (8), we can deduce
p ( y h ) p ( u h ) 0 , Ω 2 B ( p ( y h ) p ( u h ) , p ( y h ) p ( u h ) ) = ( y h y ( u h ) , p ( y h ) p ( u h ) ) 0 , Ω y h y ( u h ) 0 , Ω p ( y h ) p ( u h ) 0 , Ω .
Thus, we can get p ( y h ) p ( u h ) 0 , Ω y h y ( u h ) 0 , Ω . Using the results of [23], we  obtain
p h p ( y h ) supg h 3 / 2 , y h y ( u h ) supg h 3 / 2 .
By the triangle inequality and the relationship between · SUPG and · supg , we arrive at
p h p ( u h ) 0 , Ω p h p ( y h ) 0 , Ω + p ( y h ) p ( u h ) 0 , Ω p h p ( y h ) SUPG + y h y ( u h ) SUPG p h p ( y h ) supg + y h y ( u h ) supg h 3 / 2 .
From Theorem 1 and Equation (11), we have the following estimate:
u u h 0 , Ω h 3 / 2 + h 2 .
Inserting the above estimate into the estimates of state and adjoint state yields the final result.    □

5. 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, Figure 1b, Figure 1c in Figure 1, respectively.
The projected gradient algorithm is given below (Algorithm 1):
Algorithm 1: Projected gradient algorithm.
Require:
      Regularization parameter γ and tolerance error η .
Ensure:
  • Choose the initial value u h . Set e r r o r = 1 .
  • While e r r o r > η do
  • Solving the state equation in the discrete first order optimality system (6) to get state variable  y h ;
  • Solving the adjoint state equation in the discrete first order optimality system (6) to obtain adjoint state variable p h ;
  • Control variable u h n e w are obtained by using the projection P U a d ;
    u h n e w | E = P U a d ( 1 γ Π 1 0 ( p h | E ) ) , E T h .
  • Calculate the e r r o r :
    e r r o r = u h u h n e w L 2 .
  • Update control variable u h = u h n e w ;
  • 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].
Example 1.
Consider the optimal control Equations (1) and (2) on the unit square Ω = [ 0 , 1 ] × [ 0 , 1 ] . Let u a = 0.3 , u b = 0 , γ = 1 , and
ε = 10 10 , β = 1 3 + 10 y ( x + y 2 ) 4 1 2 5 ( x + y 2 ) 4 , δ = 1 .
The exact solutions are chosen to be
y ( x 1 , x 2 ) = 200 x 1 x 2 ( 1 x 1 ) ( 1 x 2 ) ( x 1 3 / 5 ) ( x 2 3 / 5 ) , p ( x 1 , x 2 ) = 80 x 1 x 2 ( 1 x 1 ) ( 1 x 2 ) ( x 1 3 / 5 ) ( x 2 3 / 5 ) , u ( x 1 , x 2 ) = P U a d ( p ) .
f and y d can be determined from the exact solutions y , p , u .
We choose C τ = 0.8 and let r 1 = max E T h h E β E , r 2 = min E T h h E 2 ε , r 3 = C τ δ under a fixed mesh. In order to ensure the dominance of the convection term, i.e., the choice of parameter τ E is h E β E , numerical comparisons of r 1 , r 2 , and r 3 are shown in Table 1 and Table 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).
Example 2.
We consider the Equations (1) and (2) with a modified objective functional
J ( y , u ) : = 1 2 y y d 2 + γ 2 u u 0 2
and ε = 10 4 , β = ( 1 , 0 ) , δ = 1 , γ = 1 on the unit square Ω = [ 0 , 1 ] × [ 0 , 1 ] . The exact solution of the optimal control problem is as follows:
y = 4 e ( ( ( x 1 0.7 ) 2 + ( x 2 0.7 ) 2 ) / ε ) s i n ( π x 1 ) s i n ( π x 2 ) , p = e ( ( ( x 1 0.7 ) 2 + ( x 2 0.7 ) 2 ) / ε ) s i n ( π x 1 ) s i n ( π x 2 ) , u = max { 0 , c o s ( π x 1 ) c o s ( π x 2 ) 1 } .
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 r 1 = max E T h h E β E , r 2 = min E T h h E 2 ε , r 3 = C τ δ 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.
Example 3.
In this example, we set Ω = [ 0 , 1 ] 2 , β = ( 2 , 3 ) T , δ = 1 , γ = 1 , and ε = 10 4 . The exact solutions are given by
y = 1 1 + e ( x 1 1 ) 2 + ( x 2 1 ) 2 0.7 ε , p = 1 1 + e x 1 2 + x 2 2 0.5 ε , u = max { 0.8 , p } .
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 1 2 + 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.

6. 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]).

Author Contributions

Q.W. reviewed relevant studies and the literature, completed the formula derivation and numerical simulations, and drafted the article. Z.Z. contributed the general idea of the study, checked the draft, and provided feasible suggestions and a critical revision of the manuscript. Both authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Natural Science Foundation of Shandong Province (Grants No. ZR2016JL004 and ZR2017MA020).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhu, J.; Zeng, Q.C. A mathematical formulation for optimal control of air pollution. Sci. China Ser. D. 2008, 46, 994–1002. [Google Scholar] [CrossRef]
  2. Martínez, A.; Rodríguez, C.; Vázquez-Méndez, M.E. Theoretical and numerical analysis of an optimal control problem related to wastewater treatment. SIAM J. Control Optim. 2000, 38, 1534–1553. [Google Scholar] [CrossRef]
  3. Hughes, T.J.R.; Rooks, A. Streamline upwind/Petrov Galerkin formulations for the convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Engrg. 1982, 54, 199–259. [Google Scholar]
  4. Zhou, Z.J.; Yan, N.N. A survey of numerical methods for convection-diffusion optimal control problems. J. Numer. Math. 2014, 22, 61–85. [Google Scholar] [CrossRef]
  5. Sun, T.J. Discontinuous Galerkin finite element method with interior penalties for convection diffusion optimal control problem. Int. J. Numer. Anal. Model. 2010, 7, 87–107. [Google Scholar]
  6. Leykekhman, D.; Heinkenschloss, M. Local error analysis of discontinuous Galerkin methods for advection-dominated elliptic linear-quadratic optimal control problems. SIAM J. Numer. Anal. 2012, 50, 2012–2038. [Google Scholar] [CrossRef]
  7. Xu, Q.J.; Zhou, Z.J. A mixed discontinuous Galerkin approximation of time dependent convection diffusion optimal control problem. J. Math. 2017, 2017, 6901467. [Google Scholar] [CrossRef]
  8. Fu, H.F.; Rui, H.X. A priori error estimates for optimal control problems governed by transient advection-diffusion equations. J. Sci. Comput. 2009, 38, 290–315. [Google Scholar] [CrossRef]
  9. Fu, H.F.; Rui, H.X.; Zhou, Z.J. A posteriori error estimates for optimal control problems constrained by convection-diffusion equations. Front. Math. China 2016, 11, 55–75. [Google Scholar] [CrossRef]
  10. Wang, F.Y.; Zhang, Z.Q.; Zhou, Z.J. A spectral Galerkin approximation of optimal control problem governed by fractional advection-diffusion-reaction equations. J. Comput. Appl. Math. 2021, 386, 113233. [Google Scholar] [CrossRef]
  11. Beirão da Veiga, L.; Brezzi, F.; Cangiani, A.; Manzini, G.; Marini, L.D.; Russo, A. Basic principles of virtual element methods. Math. Models Methods Appl. Sci. 2013, 23, 199–214. [Google Scholar] [CrossRef]
  12. Ahmad, B.; Alsaedi, A.; Brezzi, F.; Marini, L.D.; Russo, A. Equivalent projectors for virtual element methods. Comput. Math. Appl. 2013, 66, 376–391. [Google Scholar] [CrossRef]
  13. Beirão da Veiga, L.; Brezzi, F.; Marini, L.D.; Russo, A. Virtual element method for general second order elliptic problems on polygonal meshes. Math. Models Methods Appl. Sci. 2016, 26, 729–750. [Google Scholar] [CrossRef]
  14. Cangiani, A.; Manzini, G.; Sutton, O.J. Conforming and nonconforming virtual element methods for elliptic problems. IMA J. Numer. Anal. 2017, 37, 1317–1357. [Google Scholar] [CrossRef]
  15. Vacca, G.; Beirão da Veiga, L. Virtual element methods for parabolic problems on polygonal meshes. Numer. Methods Partial. Differ. Equ. 2015, 31, 2110–2134. [Google Scholar] [CrossRef]
  16. Antonietti, P.F.; Beirão da Veiga, L.; Mora, D.; Verani, M. A stream virtual element formulation of the Stokes problem on polygonal meshes. SIAM J. Numer. Anal. 2014, 52, 386–404. [Google Scholar] [CrossRef]
  17. Cangiani, A.; Gyrya, V.; Manzini, G. The nonconforming virtual element method for the Stokes equations. SIAM J. Numer. Anal. 2016, 54, 3411–3435. [Google Scholar] [CrossRef]
  18. Benedetto, M.; Berrone, S.; Borio, A.; Pieraccini, S.; Scialò, S. A hybrid mortar virtual element method for discrete fracture network simulations. J. Comput. Phys. 2016, 306, 148–166. [Google Scholar] [CrossRef]
  19. Berrone, S.; Borio, A. Orthogonal polynomials in badly shaped polygonal elements for the virtual element method. Finite Elem. Anal. Des. 2017, 129, 14–31. [Google Scholar] [CrossRef]
  20. Manzini, G.; Cangiani, A.; Sutton, O.J. The Conforming Virtual Element Method for the Convection-Diffusion-Reaction Equation with Variable Coeffcients; Los Alamos National Lab (LANL): Los Alamos, NM, USA, 2014. [Google Scholar] [CrossRef]
  21. Benedetto, M.F.; Berrone, S.; Borio, A.; Pieraccini, S.; Scialò, S. Order preserving SUPG stabilization for the virtual element formulation of advection-diffusion problems. Comput. Methods Appl. Mech. Eng. 2016, 293, 18–40. [Google Scholar] [CrossRef]
  22. Berrone, S.; Borio, A.; Manzini, G. SUPG stabilization for the nonconforming virtual element method for advection–diffusion-reaction equations. Comput. Methods Appl. Mech. Eng. 2018, 340, 500–529. [Google Scholar] [CrossRef]
  23. Beirão da Veiga, L.; Dassi, G.; Lovadina, C.; Vacca, G. SUPG-stabilized virtual elements for diffusion-convection problems: A robustness analysis. arXiv 2020, arXiv:2012.01104. [Google Scholar]
  24. Talischi, C.; Paulino, G.H.; Pereira, A.; Menezes, I.F.M. PolyMesher: A general-purpose mesh generator for polygonal elements written in Matlab. Struct. Multidiscip. Optim. 2012, 45, 309–328. [Google Scholar] [CrossRef]
  25. Chen, L.J.; Li, M.Z.; Xu, Q. Sinc-Galerkin method for solving the time fractional convection-diffusion equation with variable coefficients. Adv. Differ. Equ. 2020, 2020, 504. [Google Scholar] [CrossRef]
  26. Li, L.Y.; Jiang, Z.W.; Yin, Z. Compact finite-difference method for 2D time-fractional convection-diffusion equation of groundwater pollution problems. Comput. Appl. Math. 2020, 39, 142. [Google Scholar] [CrossRef]
  27. Zhang, C.Y.; Liu, H.P.; Zhou, Z.J. A priori error analysis for time-stepping discontinuous Galerkin finite element approximation of time fractional optimal control problem. J. Sci. Comput. 2019, 80, 993–1018. [Google Scholar] [CrossRef]
Figure 1. Three mesh types. (a) Hexagonal mesh; (b) Distorted square mesh; (c) Lloyd mesh.
Figure 1. Three mesh types. (a) Hexagonal mesh; (b) Distorted square mesh; (c) Lloyd mesh.
Entropy 23 00723 g001
Figure 2. The profiles of the exact solutions of state (a), adjoint state (b), control (c), SUPG-stabilized discretized optimal state y h (d), adjoint state p h (e), and control u h (f). (gi) are the profiles of the unstabilized numerically computed state, adjoint state, and control for Example 1 on Lloyd mesh with ε = 10 10 .
Figure 2. The profiles of the exact solutions of state (a), adjoint state (b), control (c), SUPG-stabilized discretized optimal state y h (d), adjoint state p h (e), and control u h (f). (gi) are the profiles of the unstabilized numerically computed state, adjoint state, and control for Example 1 on Lloyd mesh with ε = 10 10 .
Entropy 23 00723 g002
Figure 3. The convergence history of the L 2 errors, H 1 errors, standard SUPG norm errors, and the choice of parameter τ E on distorted square meshes (a), hexagonal meshes (b), and Lloyd meshes (c) with SUPG term for Example 2.
Figure 3. The convergence history of the L 2 errors, H 1 errors, standard SUPG norm errors, and the choice of parameter τ E on distorted square meshes (a), hexagonal meshes (b), and Lloyd meshes (c) with SUPG term for Example 2.
Entropy 23 00723 g003
Figure 4. The contour-lines and profiles of the SUPG-stabilized VEM discretized optimal state y h , control u h , and unstabilized state y h , control u h for Example 2 on Lloyd mesh. (a) The contour-line of stabilised state; (b) The profile of stabilised state; (c) The profile of unstabilised state; (d) The contour-line of stabilised control; (e) The profile of stabilised control; (f) The profile of unstabilised control.
Figure 4. The contour-lines and profiles of the SUPG-stabilized VEM discretized optimal state y h , control u h , and unstabilized state y h , control u h for Example 2 on Lloyd mesh. (a) The contour-line of stabilised state; (b) The profile of stabilised state; (c) The profile of unstabilised state; (d) The contour-line of stabilised control; (e) The profile of stabilised control; (f) The profile of unstabilised control.
Entropy 23 00723 g004
Figure 5. The contour-lines and profiles of the SUPG-stabilized VEM discretized optimal state y h , adjoint state p h , and control u h for Example 3 on distorted square mesh. (a) The contour-line of stabilised state; (b) The contour-line of stabilised adjoint state; (c) The contour-line of stabilised control; (d) The profile of stabilised state; (e) The profile of stabilised adjoint state; (f) The profile of stabilised control.
Figure 5. The contour-lines and profiles of the SUPG-stabilized VEM discretized optimal state y h , adjoint state p h , and control u h for Example 3 on distorted square mesh. (a) The contour-line of stabilised state; (b) The contour-line of stabilised adjoint state; (c) The contour-line of stabilised control; (d) The profile of stabilised state; (e) The profile of stabilised adjoint state; (f) The profile of stabilised control.
Entropy 23 00723 g005
Figure 6. The convergence history of the L 2 errors, H 1 errors, standard SUPG norm errors, and the choice of parameter τ E on distorted square meshes (a), hexagonal meshes (b), and Lloyd meshes (c) with the SUPG term for Example 3.
Figure 6. The convergence history of the L 2 errors, H 1 errors, standard SUPG norm errors, and the choice of parameter τ E on distorted square meshes (a), hexagonal meshes (b), and Lloyd meshes (c) with the SUPG term for Example 3.
Entropy 23 00723 g006
Table 1. Example 1: Errors and convergence rates of y and p in standard convective norm on three meshes with ε = 10 10 .
Table 1. Example 1: Errors and convergence rates of y and p in standard convective norm on three meshes with ε = 10 10 .
Meshh r 1 r 2 r 3 c y Rate c p Rate
Distorted square mesh 0.131 0.193 6.225 × 10 7 0.800 3.097 × 10 1 1.237 × 10 1
0.098 0.153 3.412 × 10 7 0.800 2.049 × 10 1 1.44 8.197 × 10 2 1.43
0.049 0.081 8.531 × 10 6 0.800 7.737 × 10 2 1.41 3.088 × 10 2 1.41
0.033 0.054 3.892 × 10 6 0.800 4.114 × 10 2 1.56 1.641 × 10 2 1.56
0.022 0.036 1.729 × 10 6 0.800 2.119 × 10 2 1.64 8.424 × 10 3 1.64
Hexagonal mesh 0.125 0.208 6.406 × 10 7 0.800 3.852 × 10 1 1.522 × 10 1
0.063 0.104 1.602 × 10 7 0.800 1.533 × 10 1 1.33 6.200 × 10 2 1.30
0.042 0.069 7.118 × 10 6 0.800 8.804 × 10 2 1.37 3.557 × 10 2 1.37
0.036 0.059 5.230 × 10 6 0.800 7.107 × 10 2 1.40 2.864 × 10 2 1.40
0.031 0.052 3.996 × 10 6 0.800 5.894 × 10 2 1.41 2.375 × 10 2 1.41
Lloyd mesh 0.201 0.274 2.136 × 10 8 0.800 7.215 × 10 1 2.954 × 10 1
0.097 0.158 5.442 × 10 7 0.800 2.539 × 10 1 1.44 1.005 × 10 1 1.49
0.070 0.116 2.603 × 10 7 0.800 1.507 × 10 1 1.55 6.047 × 10 2 1.51
0.049 0.073 1.274 × 10 7 0.800 9.224 × 10 2 1.41 3.691 × 10 2 1.42
0.035 0.052 6.429 × 10 6 0.800 5.382 × 10 2 1.66 2.145 × 10 2 1.67
Table 2. Example 1: Errors and convergence rates of y , p and u in L 2 , H 1 norms and L 2 norm, respectively, on three meshes with ε = 10 10 .
Table 2. Example 1: Errors and convergence rates of y , p and u in L 2 , H 1 norms and L 2 norm, respectively, on three meshes with ε = 10 10 .
Mesh e y , 0 Rate e y , 1 Rate e p , 1 Rate e p , 0 Rate e u , 0 Rate
Distorted square mesh 3.425 × 10 2 5.752 × 10 1 2.277 × 10 1 1.355 × 10 2 8.982 × 10 3
1.952 × 10 2 1.95 4.289 × 10 1 1.02 1.670 × 10 1 1.08 7.238 × 10 3 2.18 4.787 × 10 3 2.19
4.374 × 10 3 2.16 2.039 × 10 1 1.07 8.044 × 10 2 1.05 1.548 × 10 3 2.23 1.016 × 10 3 2.24
1.783 × 10 3 2.21 1.321 × 10 1 1.07 5.259 × 10 2 1.05 6.499 × 10 4 2.14 4.235 × 10 4 2.16
7.613 × 10 4 2.10 8.287 × 10 2 1.15 3.304 × 10 2 1.15 2.679 × 10 4 2.18 1.715 × 10 4 2.23
Hexagonal mesh 3.910 × 10 2 7.071 × 10 1 2.869 × 10 1 1.780 × 10 2 1.247 × 10 2
7.108 × 10 3 2.46 3.388 × 10 1 1.06 1.368 × 10 1 1.07 3.470 × 10 3 2.36 2.213 × 10 3 2.49
2.851 × 10 3 2.25 2.268 × 10 1 0.99 9.127 × 10 2 1.00 1.427 × 10 3 2.19 9.084 × 10 4 2.20
2.146 × 10 3 1.84 1.949 × 10 1 0.98 7.840 × 10 2 0.99 1.042 × 10 3 2.04 6.621 × 10 4 2.05
1.678 × 10 3 1.86 1.710 × 10 1 0.99 6.870 × 10 2 1.00 7.992 × 10 4 2.00 5.072 × 10 4 2.01
Lloyd mesh 9.223 × 10 2 1.084 × 10 0 4.440 × 10 1 3.715 × 10 2 2.365 × 10 2
1.853 × 10 2 2.21 4.744 × 10 1 1.14 1.900 × 10 1 1.17 7.854 × 10 3 2.14 4.781 × 10 3 2.20
8.421 × 10 3 2.35 3.303 × 10 1 1.08 1.325 × 10 1 1.07 3.661 × 10 3 2.27 2.227 × 10 3 2.28
4.151 × 10 3 2.03 2.328 × 10 1 1.01 9.237 × 10 2 1.04 1.772 × 10 3 2.09 1.139 × 10 2 1.93
2.130 × 10 3 2.06 1.729 × 10 1 0.92 6.821 × 10 2 0.93 9.518 × 10 4 1.92 6.057 × 10 4 1.95
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Back to TopTop