An HDG Method for Dirichlet Boundary Control of Convection Dominated Diffusion PDE

# An HDG Method for Dirichlet Boundary Control of Convection Dominated Diffusion PDE

Gang Chen School of Mathematics Sciences, University of Electronic Science and Technology of China, Chengdu, China (cglwdm@uestc.edu.cn).    John R. Singler Department of Mathematics and Statistics, Missouri University of Science and Technology, Rolla, MO, USA (singlerj@mst.edu).    Yangwen Zhang Department of Mathematics Science, University of Delaware, Newark, DE, USA (ywzhangf@udel.edu).
July 30, 2019
###### Abstract

We first propose a hybridizable discontinuous Galerkin (HDG) method to approximate the solution of a convection dominated Dirichlet boundary control problem. Dirichlet boundary control problems and convection dominated problems are each very challenging numerically due to solutions with low regularity and sharp layers, respectively. Although there are some numerical analysis works in the literature on diffusion dominated convection diffusion Dirichlet boundary control problems, we are not aware of any existing numerical analysis works for convection dominated boundary control problems. Moreover, the existing numerical analysis techniques for convection dominated PDEs are not directly applicable for the Dirichlet boundary control problem because of the low regularity solutions. In this work, we obtain an optimal a priori error estimate for the control under some conditions on the domain and the desired state. We also present some numerical experiments to illustrate the performance of the HDG method for convection dominated Dirichlet boundary control problems.

## 1 Introduction

Let be a Lipschitz polyhedral domain with boundary . We consider the following Dirichlet boundary control problem:

 minJ(y,u)=12∥y−yd∥2L2(Ω)+γ2∥u∥2L2(Γ),γ>0, (1.1)

subject to

 −εΔy+∇⋅(βy)+σy=fin Ω,y=uon Γ, (1.2)

where , and we make other assumptions on and for our analysis.

Researchers have performed numerical analysis of computational methods for Dirichlet boundary control problems for over a decade. Many researchers considered the standard finite element method and obtained an error estimate for the optimal control of order for all , where is the largest angle of the boundary polygon; see, e.g., [8, 37, 36]. Apel et al. in [1] considered special meshes and obtained an optimal convergence rate with . Some mixed finite element methods have also been used for Dirichlet boundary control problems because the essential Dirichlet boundary condition becomes natural, i.e., the Dirichlet boundary data directly enters the variational setting. In [22], Gong et al. used a standard mixed method to obtain an error estimate for all . Recently, we used an HDG method to obtain an optimal convergence rate for all without using higher order elements or a special mesh [30]. Moreover, the number of degrees of freedom are lower for HDG methods than standard mixed methods.

All of the above works focus on Dirichlet boundary control of the Poisson equation. However, Dirichlet boundary control problems play an important role in many applications governed by more complicated models, such as the Navier-Stokes equations; see, e.g., [25, 24, 26, 27, 21]. In order to work towards numerical analysis results for more difficult PDEs, one essential and necessary step is to fully understand the convection diffusion Dirichlet boundary control problem. Benner and Yücel in [4] used a local discontinuous Galerkin (LDG) method and they obatined an error estimate for the control of order for all . Also, very recently, we proposed a new HDG method to study this problem and obtained an optimal convergence rate for all ; see [29, 23] for more details.

However, the previous works only approximated solutions of convection diffusion Dirichlet boundary control problems in the diffusion dominated case. They did not consider the more difficult convection dominated case, i.e., . Even without the Dirichlet boundary control, solutions of convection dominated diffusion PDEs typically have layers; therefore, designing a robust numerical scheme for this problem is a major difficulty difficulty and has been considered in many works; see, e.g., [19, 32, 38, 7] and the references therein. Discontinuous Galerkin (DG) methods have proved very useful for solving convection dominated PDEs; see, e.g., [41, 13, 6, 16, 12, 33, 9] for standard DG methods and [31, 20] for HDG methods. For more information on HDG methods; see, e.g., [14, 15, 10, 11, 18, 17, 39, 40]. Moreover, there are some existing convection dominated diffusion distributed optimal control numerical analysis works; see, e.g., [3, 34, 28]. However, the techniques in the above works are not applicable for convection dominated Dirichlet boundary control problems since the solutions of (1.1)-(1.2) frequently have low regularity, i.e., with .

Formally, the optimal control and the optimal state minimizing the cost functional satisfy a mixed weak formulation of the optimality system \cref@addtoresetequationparentequation

 −εΔy+∇⋅(βy)+σy =f  in Ω, (1.3a) y =u  on Γ, (1.3b) −εΔz−∇⋅(βz)+(∇⋅β+σ)z =y−yd  in Ω, (1.3c) z =0on Γ, (1.3d) γu−ε∇z⋅n =0on Γ. (1.3e)

In this work, we use polynomials of degree to approximate the state , dual state and their fluxes and , respectively. Moreover, we also use polynomials of degree to approximate the numerical trace of the state and dual state on the edges (or faces) of the spatial mesh, which are the only globally coupled unknowns. The HDG method considered here is different from the HDG method we considered for convection diffusion Dirichlet boundary control problems in [29, 23]. A major difference is that the HDG method here has a lower computational cost.

In Section 4, we obtain an optimal convergence rate for the optimal control in 2D under certain basic assumptions on the desired state and the domain ; specifically, we prove

 ∥u−uh∥Γ≤Chs, (1.4)

for all , and the constant only depends on the exact solution, the domain and the polynomial degree. To prove the estimate (1.4), we cannot use the numerical analysis strategy from [4, 29, 23] because the constants in their error estimates may blow up as approaches zero. In order to obtain the estimate (1.4) with the constant independent of , we follow a strategy from [20] and use weighted test functions in an energy argument. However, the techniques used in [20] are not directly applicable for solutions with low regularity. Moreover, unlike all the previous Dirichlet boundary control numerical analysis works, we only assume the mesh is shape regular, not quasi-uniform. We present numerical results in Section 5 to illustrate the performance of the HDG method.

## 2 Optimality system, regularity and HDG formulation

We begin with some notation. For any bounded domain , let and denote the usual th-order Sobolev spaces on , and let , denote the norm and seminorm on these spaces. We use to denote the inner product on , and set . When , we denote , and . Also, when is the boundary of a set in , we use to replace . Bold face fonts will be used for vector Sobolev spaces along with vector-valued functions. In addition, we introduce the following space:

 H(div,Λ):={v∈[L2(Λ)]d:∇⋅v∈L2(Λ)}.

We now present the optimality system for problem (1.1)-(1.2) and give a regularity result.

### 2.1 Optimality system and regularity

Throughout the paper, we suppose is a convex polygonal domain, and let denote its largest interior angle. The optimal control is determined by the optimality system for the state and the dual state . For the HDG method, we use a mixed formulation of the optimality system; therefore we introduce the primary flux and the dual flux . The well-posedness and regularity of the mixed formulation of the optimality system is contained in the result below. The proof of Theorem 1 is omitted here since it is very similar with a proof of a similar result in [29].

###### Theorem 1.

If for some , , and the velocity vector field satisfies

 β∈[L∞(Ω)]2,∇⋅β∈L∞(Ω),σ+12∇⋅β≥0,∇∇⋅β∈[L2(Ω)]2, (2.1)

then problem (1.1)-(1.2) has a unique solution . Moreover, for any satisfying and , we have

 (u,q,p,y,z) ∈Hs(Γ)×Hs−12(Ω)×Hs+12(Ω)×Hs+12(Ω)×(Hs+32(Ω)∩H10(Ω)),

is the unique solution of \cref@addtoresetequationparentequation

 ε−1(q,r)Ω−(y,∇⋅r)Ω+⟨u,r⋅n⟩Γ =0, (2.2a) (∇⋅(q+βy),w)Ω+(σy,w)Ω =0, (2.2b) ε−1(p,r)Ω−(z,∇⋅r)Ω =0, (2.2c) (∇⋅(p−βz),w)Ω+((∇⋅β+σ)y,w)Ω =(y−yd,w)Ω, (2.2d) ⟨γu+p⋅n,v⟩Γ =0, (2.2e) for all (r,w,v)∈H(\textupdiv,Ω)×L2(Ω)×L2(Γ). Furthermore, we have Δy∈L2(Ω).

### 2.2 The HDG formulation

Let be a conforming simplex mesh that partitions the domain . For any , we let be the diameter of and denote the mesh size by . Denote the edges of by , let be the set of all edges , let be the set of edges such that , and set . Let denote the diameter of . The mesh dependent inner products are denoted by

 (w,v)Th=∑T∈Th(w,v)T,⟨ζ,ρ⟩∂Th=∑T∈Th⟨ζ,ρ⟩∂T.

We use and to denote the broken gradient and broken divergence with respect to . For an integer , denotes the set of all polynomials defined on with degree not greater than . We introduce the discontinuous finite element spaces.

 Vh :={v∈[L2(Ω)]2:v|T∈[Pk(T)]d,∀T∈Th}, Wh :={w∈L2(Ω):w|T∈Pk(T),∀T∈Th}, Moh :={μ∈L2(Eh):μ|E∈Pk(E),∀E∈Eh, and μ|Γ=0}, M∂h :={μ∈L2(E∂h):μ|E∈Pk(E),∀E∈E∂h}.

In our earlier works [29, 23], we used a local space for the spaces and . In this work, we use polynomial degree for all spaces. Since the globally coupled degrees of freedom depend on the space , the computational cost of the HDG method in this paper is much lower than the HDG method in [29, 23].

The HDG method for mixed weak form of the optimality system (2.2) is to find such that \cref@addtoresetequationparentequation

 ε−1(qh,rh)Th−(yh,∇⋅rh)Th+⟨ˆyoh,rh⋅n⟩∂Th=−⟨uh,rh⋅n⟩E∂h, (2.3a) for all rh∈Vh, (2.3b) for all (wh,ˆwoh)∈Wh×Moh, ε−1(ph,rh)Th−(zh,∇⋅rh)Th+⟨ˆzoh,rh⋅n⟩∂Th =0, (2.3c) for all rh∈Vh, (2.3d) for all (wh,ˆwoh)∈Wh×Moh, ⟨γuh+ph⋅n+τ2(zh−ˆzoh),ˆw∂h⟩ε∂h =0, (2.3e) for all ˆw∂h∈M∂h. Here, the positive stabilization functions τ1 and τ2 are chosen as τ1|∂T =∥β⋅n∥0,∞,∂T+12β⋅n+εh−1T, (2.3f) τ2|∂T =∥β⋅n∥0,∞,∂T−12β⋅n+εh−1T. (2.3g) To simplify the presentation later, we define τ=τ1+τ22=∥β⋅n∥0,∞,∂T+εh−1T. (2.3h)

### 2.3 A compact formulation

To simplify the notation, for , we denote

 B1(qh,yh,ˆyoh;rh,wh,ˆwoh) =ε−1(qh,rh)Th−(yh,∇⋅rh)Th+⟨ˆyoh,rh⋅n⟩∂Th −(wh,∇⋅qh)Th+⟨ˆwoh,qh⋅n⟩∂Th−⟨τ1(yh−ˆyoh),wh−ˆwoh⟩∂Th (2.4) B2(ph,zh,ˆzoh;rh,wh,ˆwoh) =ε−1(ph,rh)Th−(zh,∇⋅rh)Th+⟨ˆzoh,rh⋅n⟩∂Th −(wh,∇⋅ph)Th+⟨ˆwoh,ph⋅n⟩∂Th−⟨τ2(zh−ˆzoh),wh−ˆwoh⟩∂Th (2.5)

Then we can rewrite (2.3) as follows: find such that \cref@addtoresetequationparentequation

 B1(qh,yh,ˆyoh;r1,w1,ˆwo1) =−(f,w1)Th−⟨uh,τ2w1+r1⋅n⟩E∂h, (2.6a) B2(ph,zh,ˆzoh;r2,w2,ˆwo2) =−(yh−yd,w2)Th, (2.6b) ⟨ph⋅n+τ2(zh−ˆzoh),ˆw∂h⟩ε∂h =⟨γuh,ˆw∂h⟩ε∂h, (2.6c) for all (r1,w1,ˆwo1,r2,w2,ˆwo2,ˆw∂h)∈[Vh×Wh×Moh]2×M∂h.

The following basic result, which is similar to results in [29, 23], is crucial to the proof of the well-posedness of the discrete optimality system (2.3a)-(2.3e), and is also a very important part of the final stage of numerical analysis (see the proof of Lemma 14).

###### Lemma 1.

For all , we have

 B1(qh,yh,ˆyoh;rh,wh,ˆwoh)=B2(rh,wh,ˆwoh;qh,yh,ˆyoh). (2.7)
###### Proof.

Using the definitions in (2.4)-(2.5) and integration by parts give

 B1(qh,yh,ˆyoh;rh,wh,ˆwoh)−B2(rh,wh,ˆwoh;qh,yh,ˆyoh) =−⟨τ1(yh−ˆyoh),wh−ˆwoh⟩∂Th +⟨τ2(yh−ˆyoh),wh−ˆwoh⟩∂Th +(yh,β⋅∇wh)Th−⟨ˆyoh,β⋅nwh⟩∂Th−(σyh,wh)Th +(wh,β⋅∇yh)Th−⟨ˆwoh,β⋅nyh⟩∂Th+((σ+∇⋅β)wh,yh)Th =−⟨β⋅n(yh−ˆyoh),wh−ˆwoh⟩∂Th +⟨yh,β⋅nwh⟩∂Th−⟨ˆyoh,β⋅nwh⟩∂Th−⟨ˆwoh,β⋅nyh⟩∂Th =0,

where we used . This proves our result. ∎

## 3 Stability

To perform the stability and error analysis for the convection dominated boundary control problem, we need to assume some conditions on the velocity vector field and the effective reaction function .

(A1)

has a nonnegative lower bound, i.e,

 σ0:=infx∈Ω¯σ≥0. (3.1)
(A2)

has no closed curves and

 |β(x)|≠0 for all x∈Ω.
(A3)

We note that we have already assumed (A1) in Equation 2.1 in Theorem 1. We repeat the assumption here to highlight it. Also, since we are interested in the convection dominated case, (A3) is a reasonable assumption. As shown in [2], assumption (A2) implies for any integer , there exists a function such that for all , we have

 β⋅∇ψ≥2β0>0, (3.2)

where and is the diameter of . We use assumption (A3) in the analysis to remove the assumption on the meshes. Specifically, in the proofs of Lemma 11 and Lemma 15, we use assumption (A3) and a local inverse inequality to replace a global inverse inequality that has been used in all previous Dirichlet boundary control works. Therefore, we only assume is a conforming simplex partition of . All previous works on Dirichlet boundary control problems required a conforming quasi-uniform mesh. In the future, we hope to performed an a posteriori error analysis for the convection dominated boundary control problem.

###### Remark 1.

If , then assumption (A2) is the minimal known requirement that can be used to establish stability and error analysis results for numerical methods; see, e.g., [2, 20]. If instead , then we don’t need to assume (A2) and the numerical analysis is less technical. Specifically, we don’t need to prove Theorem 2 below if .

### 3.1 Preliminary material

For any nonnegative integer , we define the -projections and as follows: for any , , , , find and satisfying \cref@addtoresetequationparentequation

 (Πojv,wj)T =(v,wj)T,  ∀wj∈Pj(T), (3.3a) ⟨Π∂jq,rj⟩E =⟨q,rj⟩E,   ∀rj∈Pj(E). (3.3b)

We also define as

 ˜Π∂k|E={Π∂k|E,  E∈Eoh,0,E∈E∂h.

Then is an operator mapping to .

We first give an approximation property from [5, Theorem 4.3.8, Proposition 4.1.9], and then we prove the basic stability and approximation properties for projections.

###### Lemma 2.

Let be an integer. For any , and integer satisfying , there exists such that \cref@addtoresetequationparentequation

 |v−Im−1v|s,T ≤Chm−sT|v|m,T, (3.4a) ∥v−Im−1v∥0,∞,T ≤ChT|v|1,∞,T. (3.4b)
###### Lemma 3.

Let be a real number. For any nonnegative integer , let be a real number satisfying and let . For all , , it holds \cref@addtoresetequationparentequation

 |Πoℓv|j,T ≤C|v|j,T, ∀v∈Hj(T), (3.5a) ∥Π∂jv∥E ≤∥v∥E, ∀v∈L2(E), (3.5b) |v−Πojv|s,T ≤Chm−sT|v|m,T, ∀v∈Hm(T), 0≤s≤m, (3.5c) |v−Πo0v|ℓ,∞,T ≤Ch1−ℓT|v|1,∞,T, ∀v∈W1,∞(T), (3.5d) |v−Πojv|s,∂T ≤Chm−s−1/2T|v|m,T, ∀v∈Hm(T),0≤s+1≤m, (3.5e) ∥w∥∂T ≤Ch−1/2T∥w∥T, ∀w∈Wh. (3.5f)
###### Proof.

Equation 3.5a follows from Equation 3.5c; Equation 3.5b follows from the definition of projection; Equation 3.5e follows from Equation 3.5c and the trace inequality; and Equation 3.5f follows from the trace inequality and inverse inequality. The only thing left is to prove Equation 3.5c and Equation 3.5d.

For Equation 3.5c, in view of Equation 3.4a, an inverse inequality, and the fact that , for we have

 |v−Πojv|s,T ≤|v−Im−1v|s,T+|Im−1v−Πojv|s,T =|v−Im−1v|s,T+|Πoj(Im−1v−v)|s,T ≤|v−Im−1v|s,T+Ch−sT∥Πoj(Im−1v−v)∥0,T ≤|v−Im−1v|s,T+Ch−sT∥Im−1v−v∥0,T ≤Chm−sT∥v∥m,T.

As for Equation 3.5d, is obvious and therefore we set . By a standard scaling argument, the following stability result holds:

 ∥Πo0v∥0,∞,T≤C∥v∥0,∞,T. (3.6)

By an inverse inequality, (3.6), and (3.4b) we get

 ∥v−Πo0v∥0,∞,T ≤∥v−I0v∥0,∞,T+∥I0v−Πo0v∥0,∞,T =∥v−I0v∥0,∞,T+∥Πo0(I0v−v)∥0,∞,T ≤C∥v−I0v∥0,∞,T ≤ChT|v|1,∞,T.

In addition, we have the following super-approximation results; a similar result can be found in [2].

###### Lemma 4.

Let . Then, for any and , there holds: \cref@addtoresetequationparentequation

 ∥ηuh−Πok(ηuh)∥T ≤ChT|η|1,∞,T∥uh∥T, (3.7a) |ηuh−Πok(ηuh)|1,T ≤C|η|1,∞,T∥uh∥T, (3.7b) ∥ηuh−Πok(ηuh)∥∂T ≤Ch1/2T|η|1,∞,T∥uh∥T, (3.7c) ∥ηuh−Π∂k(ηuh)∥E ≤ChT|η|1,∞,T∥uh∥E. (3.7d)
###### Proof.

We notice that follows from , , and the trace inequality. Next, for , we have

 |ηuh−Πok(ηuh)|j,T ≤|ηuh−(Πo0η)uh|j,T+|(Πo0η)uh−Πok(ηuh)|j,T =|(η−Πo0η)uh|j,T+|Πok((Πo0η)uh−ηuh)|j,T ≤C|(η−Πo0η)uh|j,T \textupby (???), ≤C|η−Πo0η|j,∞∥uh∥T+∥η−Πo0η∥0,∞|uh|j,T ≤Ch1−jT|η|1,∞,T∥uh∥T, \textupby (???).

This proves (3.7a) and (3.7b). Similarly, for , we have

 ∥ηuh−Π∂k(ηuh)∥E ≤∥ηuh−(Πo0η)uh∥E+∥(Πo0η)uh−Π∂k(ηuh)∥E =∥(η−Πo0η)uh∥E+∥Π∂k((Πo0η)uh−ηuh)∥E ≤C∥(η−Πo0η)uh∥E \textupby (???), ≤C∥η−Πo0η∥0,∞,T∥uh∥E ≤ChT|η|1,∞,T∥uh∥E \textupby (???).

This proves (3.7d). ∎

For the analysis of the low regularity case, we need the following result from [35]:

###### Lemma 5.

If is an integer that is large enough, then there exists an interpolation operator such that for all , for all and for all , we have \cref@addtoresetequationparentequation

 (Ich(wh,ˆwoh),vh)T =(wh,vh)T, (3.8a) ⟨Ich(wh,ˆwoh),ˆvoh⟩E =⟨ˆwoh,ˆvoh⟩E, (3.8b) ∥∇Ich(wh,ˆwoh)∥Th ≤C(∥∇wh∥Th+∥h−1/2E(wh−ˆwoh)∥∂Th), (3.8c) ∥wh−Ich(wh,ˆwoh)∥Th ≤Ch(∥∇wh∥Th+∥h−1/2E(wh−ˆwoh)∥∂Th), (3.8d) where PThk+ℓ={wh∈L2(Ω):wh|T∈Pk+ℓ(T), ∀T∈Th}.