{'fraction_non_alphanumeric': 0.05765896270933474, 'fraction_numerical': 0.02565790917670815, 'mean_word_length': 3.8045977011494254, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
BDDC Algorithms for Advection-diffusion problems with HDG Discretizations

The balancing domain decomposition methods (BDDC) are originally introduced for symmetric positive definite systems and have been extended to the nonsymmetric positive definite system from the linear finite element discretization of advection-diffusion equations. In this paper, the convergence of the GMRES method is analyzed for the BDDC preconditioned linear system from advection-diffusion equations with the hybridizable discontinuous Galerkin (HDG) discretization. Compared to the finite element discretizations, several additional norms for the numerical trace have to be used and the equivalence between the bilinear forms and norms needs to be established. For large viscosity, if the subdomain size is small enough, the number of iterations is independent of the number of subdomains and depends only slightly on the sudomain problem size. The convergence deteriorates when the viscosity decreases. These results are similar to those with the finite element discretizations. Moreover, the effects of the additional primal constraints used in the BDDC algorithms are more significant with the higher degree HDG discretizations. The results of two two-dimensional examples are provided to confirm our theory.


Introduction
General hybridizable discontinuous Galerkin (HDG) methods were introduced for second order elliptic problems in [16], where the only global coupled degrees of freedom are called "numerical traces", a scalar variable. This significantly reduces the number of the degrees of freedom in the resulting global system comparing to other traditional DG methods. The HDG has been used for many different applications such as Stokes equation [17,33], Maxwell's equation [34], Helmholtz equation [35], Oseen equation [11], and advectiondiffusion equation [23,36,37]. The resulting linear system from the HDG discretization of advection-diffusion equation is nonsymmetric but usually positive definite and the generalized minimal residual methods (GMRES) can be used to solve the linear system. To reach a given accuracy, many GMRES iterations might be needed, especially for HDG discretizations with high order basis functions. Preconditioners are necessary to accelerate the convergence of GMRES. To our best knowledge, there are very few fast solvers for the advection-diffusion problems with HDG discretizations.
Domain decomposition methods are widely used preconditioner techniques for solving large sparse linear systems arising from finite element discretization of partial differential equations (PDEs). A number of domain decomposition methods have been proposed for solving advection-diffusion problems. Overlapping Schwarz methods were studied in [8][9][10], where the rate of convergence was proved to be independent of the number of subdomains if the coarse mesh is fine enough. For nonoverlapping domain decomposition methods, the standard Dirichlet and Neumann boundary conditions for the subdomain local problems cannot ensure the positive definiteness of the local bilinear forms and therefore are not appropriate, see [38,Chapter 6] and references therein for more details. A Robin boundary condition for the subdomain local problems was used in [3] to overcome this difficulty, see also [1,2]. In [41], the one-level and two-level FETI algorithms were proposed, with the same subdomain Robin boundary conditions, for solving advectiondiffusion problems.
The balancing domain decomposition by constraints (BDDC) methods, one of the most popular nonoverlapping domain decomposition methods, were introduced in [20] and analyzed in [31,32] for symmetric positive definite problems. The BDDC methods have also been extended to solving those linear systems resulting from the discretization of Stokes equations [29,48], the flow in porous media [43][44][45], and the Helmholtz equations [28]. The BDDC algorithms have also been applied to second order elliptic problems with HDG discretization [47] and the Stokes equation [48,49]. In [18], extensive numerical experiments for some additive and multiplicative BDDC algorithms with vertex constraints or edge average constraints have been studied for advection-diffusion problems with stabilized finite element discretizations. A convergence rate estimate of the GMRES was established in [46] for solving the the BDDC preconditioned advection-diffusion problems with linear finite element discretization, stabilized by the Galerkin/least squares methods [25]. The Robin boundary conditions were used for subdomain local problems and two flux-based primal constraints are added in additional to standard vertex and edge/face average constraints. Similar to [9], a perturbation approach is used to handle the asymmetry of the problem in the analysis and the partially sub-assembled finite element problem in the preconditioner is treated as a non-conforming finite element approximation.
In this paper, following the approach used in [46], we develop and analyze the BDDC algorithms for advection-diffusion problems with the HDG discretizations. Compared to the standard linear finite element discretizations, the coarse problems in the BDDC algorithms are simpler since only edge/face related constraints are needed for two/three dimensions. However, there are several difficulties in the analysis for the HDG discretizations and different technical tools are needed. First, since the BDDC preconditioned system is a reduced system for the "numerical trace" λ in the HDG discretization, we need to handle and estimate Qλ and Uλ these two operators in the HDG discretizations appropriately. Second, we need to establish the equivalent relation between the reduced norm from the bilinear form of the HDG discretization of the advection-diffusion PDE and a norm defined for λ. This norm plays a similar rule as the H 1 semi-norm in [46], where the relation is easily obtained by the definitions. Third, in the HDG discretization, since λ is defined only on the mesh element boundaries, appropriate extensions to the element interiors are needed for the error estimate of the partially sub-assembled problem in the preconditioner. A new norm, originally defined in [14], is used to replace the L 2 norm used in [46], for the non-conforming finite element approximation approach. Moreover, an error analysis for the numerical trace λ in an appropriate norm needs to be established.
The rest of the paper is organized as follows. We first describe the advection-diffusion problem and the adjoint problem in Section 2. The HDG discretizations are introduced and a reduced system for λ is formed. In Section 3, the local bilinear forms and a reduced subdomain interface problems are introduced. The BDDC preconditioner and several useful norms are provided in Sections 4 and 5, respectively. We provide an analysis of the convergence rate for our BDDC algorithms in Section 6 and all detailed proofs of the lemmas used in Section 6 are provided in Section 7. Finally, some computational results are presented in Section 8.
Here C is a positive constant independent of ǫ.

HDG discretization
Following [23,36], we introduce a new variable q = −ǫ∇u and let ρ = ǫ −1 . We obtain the following system for q and u as    ǫ −1 q = −∇u in Ω, ∇ · q + β · ∇u = f in Ω, u = 0 in ∂Ω. (2.5) Here q and u will be approximated by some discontinuous finite element spaces. We first introduce a shape-regular and quasi-uniform triangulation T h of Ω. The characteristic element size of T h is denoted by h, the element is denoted by K, and the union of faces of elements is denoted by E. The sets of the domain interior and boundary faces are denoted by E i and E ∂ , respectively.
Let Λ k = {µ ∈ M k : µ| e = 0 ∀e ∈ ∂Ω}. To make our notations simple, we drop the subscript k from now on for V k , W k , M k , and Λ k .
For each K, let (q h , u h ) ∈ (P k (K), P k (K)) such that for all K ∈ T h (ǫ −1 q h , r h ) K − (u h , ∇ · r h ) K + u h , r h · n ∂K = 0, ∀r h ∈ P k (K), (2.6) where (·, ·) K and ·, · ∂K denote L 2 -inner product for functions defined in K and ∂K, respectively. u h and q h + βu h are the numerical traces which approximate u h and −ǫ∇u h + βu h on ∂K, respectively.
The numerical trace u h = λ h for λ h ∈ Λ. ( q h + βu h ) · n takes the form: where the local stabilization parameter τ K is piecewise, nonnegative constant defined on ∂T h , see [36] for details. Following [36, Lemma 3.1], we assume where the strict inequality holds at least on one edge E in each element K.

The system for λ
Next, we reduce (2.10) to a system for λ h only by eliminating q h and u h . Defining the following operators A : The HDG method generates operator equations of the following form  In each K, given the value of λ h on ∂K, under Assumption 2.1, q h and u h can be uniquely determined, see [36,Lemma 3.1]. Namely, given λ h , the solution (q h , u h ) of (2.12) is uniquely determined by the first two equations. This elimination can be described by introducing the following operators:

By the definition, we have
AQµ Given µ ∈ M, Qµ and Uµ can be computed on each element K independently by using the value of µ on ∂K. By (2.15), we have, for all µ ∈ M and r ∈ V, (2.16) After eliminating q h and u h in each element, we could obtain a system for λ h . Once the solution λ h of (2.12) is found, the solution can be completed by computing q h and u h in each element with the given λ h . The system of λ h can be characterized by the following theorem, which is [23,Theorem 3.1]. The proof can be found in [36,Lemma 3.2].
for all µ ∈ Λ and By the definitions of the bilinear forms D, defined in (2.9), and a h , we have, a h (λ, µ) = D((Qλ, Uλ, λ) , (r, w, µ)), (2.20) for all (r, w, µ) ∈ V × W × Λ. Since a h is nonsymmetric, it is useful to write a h as the sum of the symmetric and skew-symmetric parts. In order to do this, by integration by part, we can rewrite Plug the above equation to (2.19) and we obtain For λ which is zero on ∂Ω, 1 2 β ·nλ, µ ∂T h = 0 and the symmetric and skew-symmetric parts of a h (λ, µ) are denoted by Let q, u, and λ denote the unknowns associated with q h , u h , and λ h , respectively. The matrix form of (2.12) can be written as  Eliminating q and u in each element independently from (2.23), we obtain the system for λ only We also denote A = B + Z, where B and Z are symmetric and skew-symmetric parts of A, corresponding to the bilinear forms b h and z h defined in (2.21) and (2.22), respectively. In next two sections, we will develop a BDDC algorithm to solve the system in (2.24) for the numerical trace λ.

Domain decomposition and a reduced subdomain interface problem
The domain Ω is decomposed into N nonoverlapping subdomains Ω i with diameters H i , i = 1, · · · , N, and set H = max i H i . Each subdomain is assumed to be a union of shaperegular coarse triangles and that the number of such triangles forming an individual subdomain is uniformly bounded. The subdomain interface is denoted by Γ. Γ h denotes the set of the interface nodes and can be defined as where ∂Ω i,h is the set of nodes on ∂Ω i and ∂Ω h is the set of nodes on ∂Ω. The spaces of finite element functions on Ω i are denoted by V (i) , W (i) and Λ (i) . The local bilinear forms are defined on V (i) × W (i) × Λ (i) and Λ (i) as Similarly to (2.20), by the definitions, we have, Due to the sign of the boundary integral 1 2 β · nλ (i) , µ (i) ∂T h (Ω i ) depending on β and n, the bilinear form a (i) h (·, ·) cannot be ensured to be positive definite on Λ (i) . As in [2,46], we use the Robin boundary condition on ∂Ω i and introduce The modified local bilinear forms a (i) (·, ·) are ensured to be positive definite on Λ (i) , i = 1, 2. . . . , N under the assumption (2.2). This modification will not change the global solution since λ (i) and µ (i) are continuous at element boundaries and the boundary integrals on the subdomains will cancel out between neighboring subdomains. We can write the symmetric and skew-symmetric parts of a (i) (λ (i) , µ (i) ) as In order to reduce the global problem (2.24) to a subdomain interface problem, we decompose Λ into the subdomain interior and interface parts. Let Λ Γ denote the degrees of freedom associated with Γ and Λ I be a direct sum of subdomain interior degrees of and the global problem (2.24) can be written as λ I ∈ Λ I and Λ Γ ∈ Λ Γ , We also denote the subdomain local interface numerical trace space by Λ The subdomain problems can be written as where A (i) is corresponding to the subdomain bilinear form a (i) , defined in (3.3).
We define the subdomain Schur complement S (i) We can assemble the subdomain local Schur complement S (i) Γ to obtain the global interface problem: find λ Γ ∈ Λ Γ , such that We note that S Γ is defined on the interface space Λ Γ . It is non-symmetric but positive definite. A BDDC preconditioned generalized minimal residual method (GMRES) [40] will be proposed in next section to solve (3.9).

The BDDC preconditioner
In order to introduce the BDDC preconditioner, we further decompose the subdomain interface variables to the primal and remaining variables, denoted by Λ Π and Λ ∆ , respectively. The primal variable Λ Π is spanned by subdomain interface edge/face basis functions with constant values at the nodes of the edge/face for two/three dimensions. These variables can be explicitly presented by a change of variables, see [30] and [26]. The space Λ ∆ is the direct sum of the Λ (i) ∆ , which have a zero average over each edge/face. We define Λ Γ , the partially assembled interface space, as In order to present the BDDC preconditioner, we also need some restriction, extension, and scaling operators between different spaces. The restriction operators R ∆ and R ΓΠ restricts Λ Γ to its subspace Λ Π . R Γ is the direct sum of the R (i) Γ and R Γ is the direct sum of R ΓΠ and R (i) ∆ . Recall ρ = ǫ −1 and we assume that the ρ i (x) is constant in each subdomain. Let N x be the set of indices j of the subdomains such that x ∈ ∂Ω j . A positive scaling factor δ † i (x) can be defined as follows: for γ ∈ [1/2, ∞), D,∆ is given by multiplying each row of R (i) ∆ , with the scaling factor δ † i (x) and the direct sume of R ΓΠ and R (i) D,∆ gives R D,Γ . We also denote by F Γ , the right hand side space corresponding to Λ Γ . The same restriction, extension, and scaled restriction operators will be used for the space F Γ as for Λ Γ .
With the help of R Γ , we can define the interface numerical trace Schur complement S Γ , on the partially assembled interface numerical trace space Λ Γ , as where S Γ is a direct sum of S (i) Γ , the subdomain local Schur complement defined in (3.8). The BDDC preconditioner is We have obtained the preconditioned BDDC algorithm for solving the global interface problem (3.9): find λ Γ ∈ Λ Γ , such that A GMRES method is used to solve (4.3) since the preconditioned problem is positive definite but nonsymmetric. In each iteration, one needs to multiply S Γ and S −1 Γ by a vector, which require solving subdomain Dirichlet boundary problems, Robin boundary problems, and a coarse level problem, see [30] for details. After obtaining the interface solution λ Γ , the solution can be completed by solving subdomain Dirichlet problems to find λ I .

A partial assembled finite element space, bilinear forms, and norms
In order to analyze the convergence of the GMRES algorithm for solving (4.3), in this section, we introduce several bilinear forms and corresponding norms. Similarly to [46], we first introduce a partially sub-assembled finite element space Λ, which is defined as By the definition, we can see that functions in the space Λ have a continuous primal part but in general a discontinuous dual part. Λ ⊂ Λ and the injection operator from Λ to Λ is denoted by R.
We know that the subdomain bilinear forms b (i) (·, ·), i = 1, 2, ..., N, are symmetric positive semi-definite on Λ (i) and define a semi-norm λ (i) 2 , for any w ∈ Λ, where λ (i) and w (i) are the restrictions of λ and w to Subdomain Ω i , respectively. We note that the · B and · B are norms defined for Λ and Λ, respectively.
Given λ Γ ∈ Λ Γ , we define a discrete extension of λ Γ to the interior of subdomains by By the definition, we know that λ A,Γ is obtained from λ Γ by solving subdomain advectiondiffusion problems with Dirichlet boundary conditions. We note that λ A,Γ are also well defined for λ Γ ∈ Λ Γ , and as a result λ A,Γ ∈ Λ. However, λ A,Γ does not have any energy minimization property as usual harmonic extension from symmetric elliptic problems.
Using the extension, we can define two bilinear forms for vectors in Λ Γ and Λ Γ respectively by Here p, q M represents the product q T Mp, for any given matrix M and vectors p and q. By the definitions (5.1), (5.2), and (5.3), we have the following lemma, which is similar to [46, Lemma 7.2].
Lemma 5.1. For any µ ∈ Λ, denotes its restriction to Γ by µ Γ ∈ Λ Γ . Then for any The same results also hold for functions and the corresponding bilinear forms in the space Λ Γ .
Using Lemma 5.1, for elements in the spaces Λ Γ and Λ Γ , we can define B Γ -and B Γnorms respectively as: We also introduce several useful norms, which are defined in [14,15,24]. For any domain D and λ ∈ Λ(D), we denote the L 2 norm by · D and define where |K| and |∂K| denote the n-and (n − 1)-dimensional measures of K and ∂K, respectively. Define When the domain D is Ω, we drop Ω as a subscript in notation. Namely we will use · h and · denote · h,Ω and · Ω , respectively.

Convergence rate of the GMRES iterations
In our analysis, to make the notations and proof simpler, we will focus on two dimensional case. The approach can be extended to three dimensions as well.
We need the following assumptions: Assumption 6.1. Following [4,23], we assume β ∈ W 1,∞ (Ω) satisfies: 1. β has no closed curves and β = 0, ∀x ∈ Ω; Assumption 6.2. For each face E k , the coarse level primal subspace Λ Π contains one face average, one face flux weighted average, and one face flux weighted first moment degrees of freedom such that for any w ∈ W , respectively, are the same (with a difference of factor −1 corresponding to opposite normal directions) for the two subdomains Ω i that share E k .
In this paper, we focus on the dependence of the GMRES convergence rate on the viscosity ǫ, the subdomain size H, and the mesh size h and always use c and C to represent generic positive constants independent of ǫ, H, and h. We consider 0 < ǫ ≤ 1. When ǫ is close to 1, the system (2.1) is diffusion dominant, and when ǫ is close to 0, the system (2.1) is advection dominant.
We define two constants which will be used in our analysis: as [14]. Recall C K is defined in Assumption 6.1.
We denote the preconditioned operator R T D,Γ S −1 Γ R D,Γ S Γ in (4.3) by T . The following theorem, proved in [21], gives the estimate of the convergence rate of the GMRES iteration.
Theorem 6.5. Let c and C 2 be two positive parameters such that Then where r m is the residual of the GMRES iteration at iteration m.
The following theorem is our main result for the convergence of the GMRES iterations for solving (4.3).
. C ED and C h are defined in Lemmas 6.9 and 6.11, respectively. Remark 6.7. Compare to the linear finite element discretization, studied in [46,Theorem 7.15], our upper bound has one more 1 ǫ in C ED and two additional factors ν ǫ,τ K ,h and γ h,τ K , defined in (6.1) and (6.2), respectively. For the detailed explanation, see Remark 7.21. For the common choices of τ K , γ h,τ K is usually a bounded constant. However, ν ǫ,τ K ,h can be quite large if ǫ is very small unless h is sufficiently small to balance. ν ǫ,τ K ,h appears in the lower bound as well. For the case ǫ = O(1), in Theorem 6.6, C u = 1 + log H h 2 γ h,τ and c l = 1 − CH 1 + log H h 2 , which will be positive for appropriate choice of H and h. In this case Theorem 6.5 applies and the convergence rate of our BDDC algorithm is independent of the number of subdomains and depends on H/h slightly. Due to a different proof of Lemma 6.11, our lower bound is better than that in [46,Theorem 7.15], see Remark 7.28 for details. The convergence rate deteriorates with a decrease of ǫ. For the advection-dominated cases, the performance of our BDDC algorithms are satisfactory for lower order discretization (degree 0) in our numerical experiments. The extra flux-based constraints used in our BDDC algorithms improve the performance significantly for degree 1 and 2 discretizations. However, when ǫ is very small, in order to make c l in Theorem 6.6 positive, H and h have to be tiny which might not be practical.
Similarly to the upper bound estimate for the BDDC algorithm applied to symmetric positive definite problems, for any function µ Γ ∈ Λ Γ , we define the average operator as which computes an average of µ Γ across Γ. For symmetric positive definite problems, the stability analysis of E D,Γ can be found in [26,27,47]. See [46] for nonsymmetric positive definite problems and [28,29,48] for symmetric indefinite problems.
We will provide proofs for the following Lemmas 6.8, 6.9, and 6.11 in next section.
Lemma 6.9. There exists a positive constant C such that for all µ Γ ∈ Λ Γ , The following lemma is similar to [46, Lemmas 7.9 and 7.10]. Moreover, From Lemmas 6.9, 6.8, and (6.7), we have ✷ For the lower bound, we need to estimate the L 2 type error.
Now we are ready to prove our main theorem Theorem 6.6. The proof is similar to that for [46,Theorem 7.15].
Proof of Theorem 6.6: The upper bound Equation (6.5) is given by (6.8). To prove the lower bound Equation (6.6), for any µ Γ ∈ Λ Γ , using the identity R T Γ R D,Γ = I, we have where we use Lemmas 6.8 and 6.10 for the last two steps. Canceling the common factor and using Lemma 5.1, we have where we use Lemma 6.11 for the last step. The second term in the right hand side can be combined with the left hand side and Equation (6.6) is proved. ✷

Some estimates and results
In this section, we will provide the proofs for Lemmas 6.8, 6.9, and 6.11. First, we list two useful results. For shape regular partition T h as detailed in [51], the following trace and inverse inequalities hold; see [50].
Lemma 7.1. (Trace Inequality) There exists a constant C such that where g ∈ H 1 (K), and K is an element of T h with e as an edge/face.
for any piecewise polynomial g of degree k on T h .
In order to prove Lemmas 6.8 and 6.9, we first need to establish the relation between the bilinear form a h (λ, λ) defined in (2.19) and the λ defined in (5.5) for all λ ∈ Λ. The difficulty is to estimate Qλ and Uλ these two operators. Similar estimate for the standard finite element discretization studied in [46] is much more straightforward. We will use some techniques developed in [14,15].

Estimate for the norms
We introduce two operators J Qµ and J Uµ such that for all r ∈ V and w ∈ W , We note that J Qµ and J Uµ satisfy the same equations as Qµ and Uµ defined in [14,Equation (2.9)] with τ K in [14] replaced by τ K − 1 2 β · n. By (2.16), we have We have the following results: Proof. The first and second inequalities basically follow from [14,Theorem 3.9]. However, in [14] the dependence of ǫ is not estimated explicitly. We follow [14] and give the detailed proof of these inequalities. By [14, Lemma 3.2], we have where U RT and Q RT are defined in [14, (3.5a) and (3.5b)], which are the corresponding U and Q for the elliptic problem with the RT discretization.
Following the proof of [14, Equation (3.11)], we have and therefore Following the proof of [24, Theorem 2.2], we have and the first inequality is proved.
The second inequality is obtained by using [14,Equation (3.10)]. We have The third inequality can be obtained by following the proof of [14,Equation (3.12)]. We have where we use (7.5) for the second inequality and (7.6) for the third inequality.
The fourth inequality can be proved as follows. Recall that m K (µ) is the average of µ on ∂K, defined in (5.6). We first prove that J Um Therefore J Um K (µ) = m K (µ) on K and ∇J U(m K (µ)) = 0. By the inverse inequality, the third inequality in this lemma, (5.4), and (5.5), we have Proof. See [14,Lemmas 3.3]. ✷ Lemma 7.5.
Proof. By triangle inequality we have Lemma 7.3 gives J Uµ K ≤ C µ h,K and we only need to estimate Uµ − J Uµ K . Let r = (I − J )Qµ, w = (I − J )Uµ in (7.4), then we have Applying the divergence theorem for the second term in the right hand side of the last equality, we have which implies that (7.7) By Lemma 7.4, we can obtain that Therefore, we can estimate Uµ − J Uµ K for two cases which depend on the values of h B t (Uµ − J Uµ) K and min E⊂∂K h where the second equality follows from the first equation of (7.4). By (7.7) we have where the first term of second inequality follows from trace inequality and the second term of the second inequality follows from inverse inequality, the third inequality follows from Lemma 7.3 and (7.9). Cancelling the common factor, we have Plugging in the above estimate to (7.9), we have Case II: By Assumption 6.1 and the second inequality in Lemma 7.3, we have Using the definition of B t defined in (2.11) and (7.7), we can write By (7.12) and applying Cauchy-Schwartz inequality, we have By Assumption 6.1, we have max E⊂∂K inf x∈E (τ K − 1 2 β · n) = C K > 0 and let E s K be the edge of K where inf x∈E (τ K − 1 2 β · n) attains its maximum. By (7.13), (7.8), (7.11), (7.10), and the estimation of J Uµ K in Lemma 7.3, we have Canceling the common factor (τ K − 1 2 β · n) Therefore, by Lemma 7.4 and (7.10), Combining the two cases we can obtain the result. ✷ Lemma 7.6. For µ, which is a constant on ∂K, we have Therefore for any µ ∈ Λ, where ν ǫ,τ K ,h is defined in (6.1).
Proof. If µ is a constant on ∂K, µ K = 0. By Lemma 7.3, we have and therefore ∇J Uµ = 0 in K. By (7.14) and Assumption 6.1, we have J Uµ − µ = 0 on ∂K. Therefore, for any w ∈ W (K), by integration by part, we have Plugging this identity in the second equation in (7.4) we can obtain that Take w = Uµ − J Uµ and r = (I − J )Qµ in the above equation and first equation of (7.4), we have By Assumption 6.1, we have Plugging in the first formula of (7.4), for any r ∈ V we have which implies that ∇(I − J )Uµ| K = 0. Therefore, (I − J )Uµ is a constant on K. Due to the fact that (I − J )Uµ = 0 on ∂K, we obtain (I − J )U(µ)| K = 0. Therefore, for any µ ∈ Λ, where m K (µ) is the average of µ on ∂K defined in (5.6).
Proof. By Lemma 7.6 and [14, Lemma 3.5 (i)], we have, on ∂K, Therefore, where we use (7.16) for the second equality, the trace theorem for the third inequality, Lemma 7.5 for the fourth inequality, and (5.5) and (5.4) for the last inequality. ✷ Lemma 7.8.
Proof. Let r = (I − J )Qµ and w = (I − J )Uµ in (7.4), and we have We will estimate these five terms separately as follows: For the first term, we have where we use Lemmas 7.5 and 7.3 for the last inequality. For the second term, by Lemmas 7.5 and 7.3 we have For the fourth term, we first use integration by part for the second equality and Lemma 7.3 . For the last term, by Assumption 6.1 and Lemma 7.3, we have Combing the estimation of these five terms and we obtain Cancelling the common terms in both sides, we have . For all µ with zero boundary on quasiuniform meshes, we have Proof. See [14,Lemma 3.7]. ✷ Lemma 7.10.
Proof. By (2.21) and (2.22), we have Here we use Lemma 7.3 for the first term, Lemma 7.8 for the second and third terms, Lemma 7.3 for the fourth term, and Lemma 7.5 for the last term.
Combining the above estimate with Lemma 7.9, we obtain the upper bound estimate.
On the other hand, by the similar proof in [14, Lemma 3.8], we can obtain that ǫ −1 Qµ K ≥ C µ K , therefore ✷ In the following, we establish the relation between · h,Ω i and · Ω i for the functions in Λ (i) , which have zero edge average on ∂Ω i . We also prove similar result as Lemma 7.9 for the functions in Λ.
Let Λ 0,(i) be the zero-order numerical trace space in Ω i and Q 0 be the L 2 orthogonal projection from Λ (i) into Λ 0,(i) .
By the Poincare inequality (see [7]) and the scaling argument for P 1 non-conforming element, we have Combing this with (7.19), we have Therefore, combing (7.20) and (7.21), we have ✷ Using Lemma 7.12 and the proof of Lemma 7.10, we can prove the following equivalence for subdomain functions.
Lemma 7.13. If µ (i) ∈ Λ (i) and µ (i) is edge average zero on ∂Ω i , then We also note that the lower bound is always hold for any µ (i) ∈ Λ (i) . Assumption 6.3 assume that the subdomains form a shape regular coarse mesh of Ω, denoted by T H . We also denote by V 0 H , the P 1 non-conforming finite element space on this coarse mesh formed by the subdomains. Let Λ 0 H be the zero-order numerical trace space in Ω on this mesh. Similarly to X h defined in (7.18), we can establish an isomorphism where λ ∈ Λ 0 H and x E denotes the midpoint of element edge E of the coarse mesh. Let Q V H denotes the L 2 orthogonal projection from V 0 h to V 0 H and I V h is an intergrid transfer operator from V 0 H to V 0 h , defined in [24]. A similar operator can be found in [6]. By the definition/construction, we have that Q V H and I V h both preserve the edge average over the edges in T H .
By [24, Lemma 4.5 and Lemma 4.6], we have and We have the similar result as Lemma 7.9 for the functions in Λ.

Next we give an estimation of
Since µ ∈ Λ, the average of µ (i) are the same across the subdomain interfaces. Therefore, µ (i) 0 and X h µ (i) 0 have the same average across the subdomain interfaces, see the proof in Lemma 7.12. Given an edge E ij shared by Subdomains Ω i and Ω i , we have Due to the linearity of u H . Since µ has zero boundary condition, u H is zero on ∂Ω. By the Poincare inequality (see [7]) for P 1 non-conforming finite element functions, we have H have zero edge averages on ∂Ω i . We also have, by the Poincare inequality in [7],

Estimates for the bilinear forms
Recall that the symmetric and skew-symmetric parts of the subdomain local bilinear form a (i) (η (i) , µ (i) ), defined (3.3), are denoted by b (i) (η (i) , µ (i) ) and z (i) (η (i) , µ (i) ), defined in (3.4) and (3.5), respectively. We will give the estimates of these bilinear forms which are useful for our lower and upper bound estimates.
Lemma 7.15. There exists a positive C such that for all where we use the definition of · B (i) and the Cauchy-Schwarz inequality. We can prove the bound for z (i) (η (i) , µ (i) ) as follows: where Assumption 6.1 is used in the last step. By Lemmas 7.6 and 7.5, we have

By Lemma 7.7, we have
Plugging in all these estimates in (7.27) and use Lemma 7.10, we have Proof. By integration by part, we have Since where we use Lemmas 7.5 and 7.6 for the first term, Lemma 7.5 and (3.4) for the second term, Assumption 6.1, Lemma 7.7 and (3.4) for the third term, and Lemmas 7.1, 7.2, 7.5, Assumption 6.1, and Lemma 7.7 for the last term. We conclude the proof by using the lower bound in Lemma 7.10 for Subdomain Ω i . ✷ Remark 7.17. Compare to [46], we cannot control L 2 -norm of Uµ when the velocity field β is divergence free. Therefore, we have additional 1 √ ǫ in our estimate. Lemma 7.18. There exists a positive C such that for all η, µ ∈ Λ Proof. Using Lemma 7.15, we have where we use Lemmas 7.9 and 7.10 for the last two inequalities, respectively. By Lemma 7.15, we have Similarly, by using Lemma 7.14, we can prove for η, µ ∈ Λ.
Similarly to the proof of Lemma 7.16, by integration by part, we have where we use β · nη, µ ∂T h = 0 for µ = 0 on ∂Ω for the third equality. For the last inequality, we use Lemma 7.6 and 7.5 for the first term, Lemmas 7.5 and 7.9 for the second term, Assumption 6.1, Lemmas 7.7 and 7.1, 7.2 for the last two terms. ✷ Proof of Lemma 6.8: where we use Lemma 7.18 for the second inequality and Lemma 5.1 again for the last equality. The result for λ Γ , µ Γ ∈ Λ Γ with the corresponding norm can be proved similarly using Lemma 7.18. Similarly, by the definition (5.2), we have where we use Lemma 7.18 for the second inequality, Lemma 7.10 for the third inequality, and Lemma 5.1 for the last equality.

Estimate of the average operator
We define the following elliptic problem: −ǫ∆u = f, in Ω, u = 0, on ∂Ω.
Let a e h (η, µ) be the bilinear form defined in [47, (2.7)] and define the norm Given µ H,Γ is defined as A e(i) (7.30) and define S e(i) Γ norm as and by Lemma 7.13, we have v (i) Therefore, where we use (7.32) for the last step. Let µ where we use Lemmas 7.19 again. By Lemma 7.13, we have

✷
Remark 7.21. Due to the HDG discretization, we have the same bounds for both two and three dimensions. Compare to the result for the linear conforming finite element discretization in [46,Lemma 7.8], our bound has one more 1 ǫ and two additional factors ν ǫ,τ K ,h and γ h,τ K , defined in (6.1) and (6.2), respectively. The additional 1 ǫ is due to the estimate in Lemma 7.16, see Remark 7.17 for more details. The first factor ν ǫ,τ K ,h is from Lemma 7.6, where ∇Uµ K ≤ ν ǫ,τ K ,h µ K . However, for the linear conforming finite element, ∇µ K = |µ| H 1 (K) is hold simply by the definition. The second factor γ h,τ K is from Lemma 7.20, which is the upper bound for the average operator of the corresponding elliptic problem with the HDG discretization. For the common choices of τ K , γ h,τ K is usually a bounded constant. However, ν ǫ,τ K ,h can be quite large if ǫ is very small unless h is sufficiently small to balance.

Estimate of the L 2 type error
This subsection is devoting to proving Lemma 6.11. The usual techniques to estimate the L 2 error is a duality argument [5]. [46,Lemma 7.12], which is a standard finite element version of Lemma 6.11, was proved using these techniques. However, there are several difficulties for the HDG discretization: (1) the preconditioned system (4.3) is a reduced system for the variable λ ∈ Λ. We need to handle Qλ and Uλ these two operators appropriately; (2) λ is defined only on mesh element boundaries, we need appropriate extensions to the element interiors.
In order to overcome these difficulties, we first introduce a local extension operator S K i as that defined in [14,Section 5.2]. Given any λ ∈ Λ, let λ K be the restriction of λ on ∂K, where K is an element in T h . Let E i denote the edge of K opposite to the ith vertex of K. We define S K i λ K in P k+1 (K) by for i = 1, · · · , n + 1 all n + 1 edges of K. Define Lemma 7.22. Equations (7.33) and (7.34) uniquely define an S K i λ K in P k+1 (K). For all λ ∈ Λ, we have

38)
where Q k is a L 2 projection to P k (K) and ν ǫ,τ K ,h is defined in (6.1).
Proof. This can be proved by a modification on the proof of [14, Lemma 5.7] with Lemma 7.5 and Lemma 7.6. ✷ In our proof, we will also need to use another extension X h , defined in (7.18), to extend λ defined on the element boundary to the function defined in the element interior by connecting with the P 1 non-conforming finite element functions. Now we are ready to use the duality argument to prove Lemma 6.11.
The following lemma gives an error estimate for the HDG solution of the dual problem. Since we need to obtain the error estimate in B-norm, the results in [23] cannot be used directly.
Cancelling the common factor both sides, we obtain Using the properties of the projections, see [23], we have ✷ We will prove the following lemma which is similar to [46,Lemma 5.1]. However, there are several technical differences. Since for any q ∈ Λ, q is only defined on the element boundary, we need to use Uq which is the extension of q into interior element to define the L 2 inner product. Due to the HDG discretization, we use ϕ g , the solution of HDG, in the definition of L h , which requires an error bound provided in Lemma 7.23. This new definition of L h simplifies the proof of Lemma 7.25. Moreover, the estimation of the term related to q requires more tools. Lemma 7.24. Let Assumption 6.2 hold and L h (q, ϕ g ) = a(q, ϕ g ) − (Uq, g) T h , for q ∈ Λ. There exists a constant C, which is independent of ǫ, H, and h, such that for all q ∈ Λ, Proof: For any q ∈ Λ, we have where we use the definition of a (i) defined in (3.3) for the second equality and (3.2) for the third equality. The Galerkin-orthogonality condition for the dual discretization and the definition of D (i) (3.1) are used for the fifth equality. The definitions of Q * g , u * g , ϕ * g and (7.42) for the second, third and fourth terms in the last equality. Therefore we have L h (q, ϕ g ) = a(q, ϕ g ) − (Uq, g) Ω i (7.43) where we use ϕ * g = u * g on the element boundary for the second equality. Denote the common average of q on the edge E ij of ∂Ω i by q E ij .
where we have also subtracted the constant average values q E ij from q (i) , which does not change the sum. Then, from Assumption 6.2, we know that where I (i) H u g represents the interpolation of u g into the finite element space on the coarse subdomain mesh. Plug (7.44) in (7.43), we have We will show in the following that each of the three terms in (7.45) can be bounded by CC L u * g H 2 q B , where C is a positive constant independent of H and h. For the first term I 1 , from the Cauchy-Schwarz inequality, we have For the second factor, we can estimate as follows. Let q (i) 0 = Q 0 q (i) and we have (7.48) The first term in (7.48) can be estimated by Lemma 7.11 as The second term in (7.48) can be estimated by the isomorphism X h defined by (7.18). We follow [46,Equation (23)] and have where we have used the definition of X h for the first inequality, a trace theorem and a Poincaré inequality of P 1 non-conforming finite element functions for the second inequality, (7.19) for the third inequality, and Lemma 7.11 for the last inequality. Therefore, we have Combining (7.46), (7.47), and (7.49), we have the following bound for I 1 , where we use the Cauchy-Schwarz inequality and Lemma 7.10 in the last step.
To derive a bound for I 2 , we find from the Cauchy-Schwarz inequality that Using a trace theorem and Lemma 6.4, we have, for the first factor on the right hand side of (7.50), For I 3 , we have where the second inequality follows from (7.49).
where we have used (7.40) for the third equality.
Dividing by ϕ λ − ϕ λ B on both sides and denoting ϕ λ − ϕ λ by q, we have, by the fact that Proof: For any v ∈ R( Λ), denote its continuous interface part by v Γ ∈ R Γ ( Λ Γ ). Given µ Γ ∈ Λ Γ , from Lemma 5.1 and the fact that The following lemma is similar to [46,Lemma 7.12]. However, we have to use S-norm defined in (7.35) instead of L 2 norm, which requires more technical tools.
Proof: Given any q (i) , λ (i) ∈ Λ (i) , by (7.35), (7.34), and (7.41), we have Using the Friedrichs estimate and (7.38), we have . Combining these two estimates, we have We have, from (7.57), (7.39) and (7.40), that Since ϕ λ ∈ Λ, by Lemma 7.26, we know that a(w A,Γ − µ A,Γ , ϕ λ ) = 0. Therefore, and, by Lemmas 7.18 and 7.25, (7.58), Therefore, we have where we use Lemma 6.10 for the last step. Proof of Lemma 6.11: Proof. Given µ Γ ∈ Λ Γ , let w Γ = S −1 Γ R D,Γ S Γ µ Γ . Recall the restriction operators R Γ . To make the notation simpler, let For the first term, since w Γ ∈ Λ Γ , R T D,Γ w Γ and w Γ have the same edge average on ∂Ω i . We use Lemmas 7.12, 6.9, 6.10 and have For the second term, we use Lemmas 7.22, 7.27 and have  [46,Lemma 7.14], we estimate the first and second term separately and therefore our estimate does not have the factor H h . However, the additional factors ν 6 ǫ,τ K ,h and γ h,τ K are still in our estimate which are from the HDG discretization. An additional 1 ǫ is due to the estimate of the average operator E D .

Numerical experiments
We present some numerical results for two examples, which were considered in [41,46]. In [46], the performance of the BDDC algorithms for Test Problem I and II are very similar. Therefore, we only present the results for Test Problem I (thermal boundary layer) and III (Rotating flow field) in [46]. We decompose our domain Ω = [−1, 1]×[−1, 1] into subdomains and each subdomain into triangles. We then discretize (2.4) with HDG discretization of degree 0, 1, and 2 polynomials. In (2.7), we take the stabilizer β E = max(sup x∈E (β · n, 0)) for any E ∈ ∂K and ∀K ∈ T h , as defined in [23, (2.10)]. We take f = 0 in all the examples as in [46].
We use the GMRES method without restart and work with the L 2 -norm. The iteration is stopped when L 2 -norm of the initial residual is reduced by a factor 10 −10 . The convergence rates in K-norm are quite similar with using the L 2 -norm [46]. We present the results of our BDDC algorithms employing different primal continuity constraints. For BDDC 1 , the subdomain edge average constraints are chosen as coarse level primal variables, while in BDDC 2 , besides the subdomain edge average constraints, one additional edge flux weighted average constraint, defined in Assumption 6.2, for each edge, are also added to the coarse level subspace. For BDDC 3 , the first moment defined in Assumption 6.2, is added additional to those used in BDDC 2 .
We note that the one-level and two-level Robin-Robin methods, which were developed and extended in [1][2][3], are closely related to our BDDC algorithms. In [46], the performance of these Robin-Robin algorithms are compared with similar BDDC algorithms as we discussed in this paper. The BDDC with extra flux weighted edge constraints overperformed. We do not present any results for the these algorithms here. See [46] for more details.
In the first set of our experiments, we change the number of the subdomains with a fixed subdomain local problem. The results are reported in Table 1. We note that, for large values of the viscosity ǫ (ǫ ≥ 0.01) with different degrees of the basis functions in HDG, the iteration counts for all BDDC algorithms are independent of the number of the subdomains. However, with the increasing of the degrees of the basis functions, the system become more ill-conditioned and those additional constraints make the BDDC perform better.
With a decreasing value of ǫ (ǫ < 0.01), we can see the scalability of the algorithms deteriorates. But the more constraints used in the BDDC algorithms, the performance are better. We also note that the numbers of the iterations are bounded for all algorithms even when the viscosity ǫ goes to zero. These observations are similar to those for the FETI algorithms in [41]. When ǫ ≤ 10 −5 , the iterations are double as we half H, the size of the subdomains, for all BDDC algorithms. In our theory, we require H small enough to obtain the scalable results. In order to confirm that, we did larger experiments for those degree zero basis functions, which we can afford. The results are listed in Table 2. When we increase the number of the subdomains to 128 for ǫ = 10 −3 (ǫ ≥ 10 −4 for BDDC 3 ), the iterations do not increase anymore. However, for smaller ǫ, 128 is still not larger enough to see the scalability of the algorithms.
In the second set of our experiments, we change the size of the subdomain local problems with a fixed number of the subdomains. The results are listed in Table 3. From the results, we can see that the number of iterations increases slightly with increasing H h , the size of the subdomain local problems, when the viscosity ǫ is large (ǫ ≥ 10 −4 ). These results are consistent with the results of the domain decomposition algorithms for symmetric, positive definite problems. Moreover, the more constraints are used in the BDDC algorithms, the better performance we obtained.
However, with the decreasing of the values of ǫ, the iterations for all BDDC algorithms are almost independent of the values of ǫ, H h , and the constraints of the BDDC algorithms. These results are consistent with the results in [46].

Rotating flow field (Test Problem II)
The second test case is more difficult one, see [41,46]. The velocity field be β = (y, −x) and the boundary condition is: We also have two sets of the numerical experiments for this example. In the first set of our experiments, we change the number of subdomains with a fixed size of the subdomain local problems. The results are reported in Table 4.
Again, we note that, for large values of the viscosity ǫ ≥ 10 −3 , the iteration counts are independent of the number of the subdomains for all three BDDC algorithms with different degrees of the basis functions. However, the additional constraints are more significant in this difficult case. For the basis functions with degree 0 and 1, the iteration counts with BDDC 3 are almost independent of the number of subdomains, even with very small ǫ = 10 −6 , while this is not the case for the BDDC 1 algorithm with ǫ < 10 −4 for degree 0 basis function. The performance of BDDC 2 is between these two algorithms.
With the degree 2 basis functions, the iteration counts increase rapidly with the decreasing of ǫ. For ǫ = 10 −6 , the number of iterations for the BDDC 1 is almost 20 times of those for BDDC 3 for large number of subdomains. Even though the iteration counts for BDDC 3 is increasing in this case, it is much less than double when we half H, the size of the subdomains.
In the second set of our experiments, we change the size of the subdomain local problems with a fixed number of the subdomains. The results are reported in Table 5. From the results, we can see that the number of iterations increases with increasing H h for all values of the viscosity ǫ for all BDDC algorithms. The problems become more difficult when we increase the degrees of the HDG basis functions. However, the BDDC 3 algorithm works significantly better than the others.