A conservative hybrid method for Darcy flow

We present a hybrid mimetic spectral element formulation for Darcy flow. The discrete representations for 1) conservation of mass, and 2) inter-element continuity, are topological relations that lead to sparse matrix systems. These constraints are independent of the element size and shape, and thus invariant under mesh transformations. The resultant algebraic system is extremely sparse even for high degree polynomial basis. Furthermore, the system can be efficiently assembled and solved for each element separately.


Introduction
Hybrid formulations [1,3,9] are classical domain decomposition methods which reduce the problem of solving one global system to many small local systems. The local systems can then be efficiently solved independently of each other in parallel.
In this work we present a hybrid mimetic spectral element formulation to solve Darcy flow. We follow [8] which render the constraints on divergence of mass flux, the pressure gradient and the inter-element continuity metric free. The resulting system is extremely sparse and shows a reduced growth in condition number as compared to non-hybrid system. This document is structured as follows: In Section 2 we define the weak formulation for Darcy flow. The basis functions are introduced in Section 3. The evaluation of weighted inner product and duality pairings are discussed in Section 4. In Section 5 we discuss the formulation of discrete algebraic system. In Section 6 we present results for a test case taken from [7].

Darcy flow formulation
For Ω ∈ R d , where d is the dimension of the domain, the governing equations for Darcy flow, are given by, where, u is the velocity, p is the pressure, f the prescribed source term, A is a d × d symmetric positive definite matrix,p andû n are the prescribed pressure and flux boundary conditions, respectively. Notations For f , g ∈ L 2 (Ω), ( f , g) Ω denotes the usual L 2 -inner product.
For vector-valued functions in L 2 we define the weighted inner product by, where (· , ·) denotes the pointwise inner product. Duality pairing, denoted by ·, · Ω , is the outcome of a linear functional on L 2 (Ω) acting on elements from L 2 (Ω).
Let Ω K be a disjoint partitioning of Ω with total number of elements K, and K i is any element in Ω K , such that, K i ∈ Ω K . We define the following broken Sobolev spaces [2], H (div;

Weak formulation
The Lagrange functional for Darcy flow is defined as, The variational problem is then given by: For given f ∈ L 2 (Ω K ), find u ∈ H(div; Ω K ), p ∈ L 2 (Ω K ), λ ∈ H 1 2 (∂Ω K ), such that,

Basis functions
Primal and dual nodal degrees of freedom Let ξ j , j = 0, 1, ..., N, be the N + 1 Gauss-Lobatto-Legendre (GLL) points in I ∈ [−1, 1]. The Lagrange polynomials h i (ξ) through ξ j , of degree N, given by, form the 1D primal nodal polynomials which satisfy, h i (ξ j ) = δ i j . Let a h and b h be two polynomials expanded in terms of h i (ξ). The L 2 -inner product is then given by, are the nodal degrees of freedom. We define the algebraic dual degrees of freedom, a, such that the duality pairing is simply the vector dot product between primal and dual degrees of freedom, Thus, the dual degrees of freedom are linear functionals of primal degrees of freedom.

Primal and dual edge degrees of freedom
The edge polynomials, for the N edges between N + 1 GLL points, of polynomial degree N − 1, are defined as [4], Let p h and q h be two polynomials expanded in edge basis functions. The inner product in L 2 space is given by, are the edge degrees of freedom. As before, we define the dual degrees of freedom such that, A similar construction can be used for dual degrees of freedom in higher dimension. For construction of the dual degrees of freedom in 2D see [8] and for 3D see [10].

Differentiation of nodal polynomial representation
Let a h (ξ) be expanded in Lagrange polynomial, then Therefore, taking the derivative of a polynomial involves two steps : First, take the difference of degrees of freedom; and second, change of basis from nodal to edge [4].

Discrete inner product and duality pairing
For 2D domains, the higher dimensional primal basis are constructed using the tensor product of the 1D basis.
For the weak formulation in (2) we expand the velocity u h in primal edge basis as, (4) Weighted inner product Using (1) and the expansions in (4), the weighted inner product is evaluated as, where, u K i are the degrees of freedom in element K i , and For mapping of elements please refer to [6].

Divergence of velocity
Divergence of velocity, ∇ · u h , is evaluated using (3), but now for 2D, The pressure is expanded in the dual basis e i (ξ) e j (η). These basis are dual to the basis in which ∇ · u h is expanded in (5). Therefore the weak constraint on divergence of velocity is a duality pairing evaluated as, where E 2,1 represents the discrete divergence operator. It is an incidence matrix that is metric-free and topological, and remains the same for each element in Ω K . For an extensive discussion on the incidence matrix, see [6]. For an element of degree N = 3,

Connectivity matrix
The connectivty matrix ensures continuity of the velocity across the elements. λ is the interface variable between the elements that acts as Lagrange multiplier that imposes the constraint given by, where N is the discrete trace operator. It is a sparse matrix that consists of 1, −1 and 0 only. For construction of N please refer to [5]. E N is the assembled N for all the elements. For discretization, K = 2 × 2, N = 2, E N is shown in (6). The matrix size of E N is 8 × 64, but it has only 16 non-zero entities. It is an extremely sparse matrix that is metric-free and the location of +/-1 valued entries depend only on the connection between different elements.

Discrete formulation
Using the weighted inner product and duality pairings discussed in Section 4, we can write the discrete form of weak formulation in (2) as, where, A is an invertible block diagonal matrix given by, E N is as given in (6), is the only dense matrix and also the only component that changes with each local element, K i . E N is a sparse incidence matrix for the global system and E 2,1 is a sparse incidence matrix for the local systems that remain the same for each element.
Using the Schur complement method, the global system (7) can be reduced to solve for λ, [1], To evaluate λ in (9) we need A −1 that can be calculated efficiently by taking inverse of each block of A separately. This part can be easily parallelized. Once the λ is determined the solution in each element, K i , can be evaluated by solving for, Here the inverse of the local block in the RHS is already evaluated during (9). As the local systems are independent of each other (10) can also be evaluated separately for each element and easily parallelized. The system (9) solves for interface degrees of freedom between the elements and will always be smaller than the full global system. For a comparison of the size of λ system with full system see Table 1 for 2D systems, and Table 2 for 3D systems. In Table 1 & 2 (left) we see that, for constant K, increasing the order of polynomial basis the growth in size of λ system is less than the growth in size of full system. Thus, hybrid formulations are beneficial for high order methods, in 2D and in 3D, where local degrees of freedom of an element are much higher than interface degrees of freedom.
In Table 1 & 2 (right) we see that, for constant N, the λ system is certainly smaller than the full system, although the growth rate in size of λ and full systems does not change significantly.

Results
In this section we present the results for a test problem from [7] by solving system (7). The domain of test problem is, Ω ∈ [0, 1] 2 . The source term is defined as, , and Dirichlet boundary conditions are imposed along the entire boundary, Γ D = Γ and Γ N = ∅. We solve this problem on an orthogonal and a highly curved mesh, see  compared to the continuous element formulation, 117504. Here, the sparsity is due to use of algebraic dual degrees of freedom and is not because of hybridization of the scheme. In Fig. 3, on the left we compare the growth in condition number, for the λ only system with continuous element system, for N = 7 on the curved mesh, with increasing number of elements, K. We observe similar growth rates for hybrid and continuous formulation, however the condition number for continuous elements formulation is almost O 10 2 higher. On the right we see the growth in condition number with increasing polynomial degree for K = 9 × 9 on the curved mesh. A suppressed growth rate in condition number for hybrid formulation is observed. Thus hybrid formulations are beneficial for high order methods. In Fig. 4 we show the L 2 -error for ∇ · u h − f h . On the left side as a function of element size, h = 1/ √ K, and on the right side as a function of polynomial degree of the basis functions. In both cases the maximum error observed is of O 10 −12 .
In Fig. 5, on the top two figures we show the error in the H (div; Ω) norm for the velocity; and at the bottom two figures we show the error in L 2 (Ω) norm for the pressure. On the left we have h-convergence plots, and on the right we have N-convergence plots. In all the figures, for the same number of elements, K, and polynomial degree, N, the error is higher for the curved mesh.
On the left we see that the error decreases with the element size. The slope of error rate of convergence is N, which is optimal for both curved and orthogonal meshes. On the right we see exponential convergence of the error with increasing polynomial degree of basis for both orthogonal and curved meshes.