A High-Order HDG Method with Dubiner Basis for Elliptic Flow Problems

We propose a standard hybridizable discontinuous Galerkin (HDG) method to solve a classic problem in fluids mechanics: Darcy’s law. This model describes the behavior of a fluid trough a porous medium and it is relevant to the flow patterns on the macro scale. Here we present the theoretical results of existence and uniqueness of the weak and discontinuous solution of the second order elliptic equation, as well as the predicted convergence order of the HDG method. We highlight the use and implementation of Dubiner polynomial basis functions that allow us to develop a general and efficient high order numerical approximation. We also show some numerical examples that validate the theoretical results.


Introduction
The study of flow processes on porous media has gained particular strength during last decades. This area is relevant from the theoretical point of view but even more when it is related to realistic applications, e.g. geological and environmental sciences [1], [2], [3]. In such studies the construction of accurate mathematical models leads to a better understanding of complex systems [4]. Further, the design of accurate and strong mathematical numerical methods provide a suitable tool in the cases when experimental result are not feasible.
We consider Darcy's law for single phase flow over porous media. This empirical model relates the pressure and the velocity of an incompressible fluid through a porous medium. We remark that Darcy's can be derived from the Navier-Stokes equations via homogenization and it is valid for slow, viscous flow (see [5], [3]).
Moreover, one of the current challenges on simulating fluid flow trough porous media is related with the highly heterogeneity of the media, e.g. heterogeneous reservoirs. Therefore, we focus on the numerical solution of the heterogeneous and anisotropic case of Darcy's law, as mentioned in [6], [7]. This is a general framework from where we aim to extend classical mathematical results and the application of novel numerical methods.
Let Ω ⊆ R 2 be a bounded and simply connected domain with polygonal boundary Γ. For a fluid velocity u and preasure p, we consider the second-order elliptic mixed problem div u = f, in Ω, (1b) where K := K(x) is a symmetric and positive definite tensor that represents the permeability of the media and n is the unit outward normal to Γ. The given data include the flux source f ∈ L 2 (Ω) and the flow at the boundary g ∈ H 1/2 (Γ). The problem (1) is the mixed formulation of Darcy's law and describes the single-phase flow through a porous media. The existence of the |34 Ingeniería y Ciencia solution (u, p) ∈ H(div; Ω) × H 1 m (Ω) of (1) is provided by the compatibility condition and the uniqueness is guaranteed in the space H 1 m (Ω) := q ∈ H 1 (Ω) : We propose the application of the Hybridizable Discontinuous Galerkin (HDG) method to approximate the solution of (1). We refer to [8] for a general overview of mixed methods applied to Darcy's law. This area has been widely studied from all fronts; many works during the last 20 years have shown progress in the numerical approximation of the Darcy equation and other equations of fluid mechanics using mainly the finite volume method and finite differences. The discontinuous Galerkin methods (DG) have the advantage (when compared with classical finite element methods) of allowing handling heterogeneous media and discontinuous solutions. In [9] a discontinuous Galerkin method is proposed to solve the Navier-Stokes problem using Raviart-Thomas elements and constitutes one of the first approaches used to solve this equation based on its mixed formulation. Also, in [10], the DG formulations to solve classic equations of fluid dynamics and heat transfer are explained. An important drawback of classical DG methods is the computational cost that involves the increasing number of degrees of freedom.
The hybridizable discontinuous Galerkin methods are a family of DG methods that allows to find the approximate solution by solving an equivalent system of equations associated to the skeleton of the partition of the domain. We refer to [16] and [17] for details about the HDG formulation. With respect to fluid base problems, in [18] a detailed study of the HDG method is showed. They applied the HDG method to problems associated with compressible fluids using Runge-Kutta type methods as well as numerical differentiation for time-dependent problems. More recently, we highlight some relevant results to solve the Stokes equation and other fluid problems using HDG, in [19], [20], [21], [22]. However, many of these papers lack a complete error analysis or only consider the case of constant permeability. In general, the mathematical foundation of the HDG method, for this kind of problems, has not been as studied as other DG-type methods. These motivate us to work in this direction and it is the main contribution of this paper.
In this work we also want to exploit an advantage of the discontinuous Galerkin method: the capability to obtain high order approximations. To these end we present and implement the Dubiner polynomial basis, that allows us to work with high order polynomials without additional processing of the scheme. The Dubiner basis have been studied in [23] and [24].
The paper is organized as follows. In Section 2 we give some of the main ideas of the HDG method. In Section 3 we show a complete error analysis of the HDG method applied to the problem (1). Finally, in Section 4 the Dubiner basis are briefly commented and in Section 5 we discuss two classical numerical tests that reflects the super-convergence of the HDG scheme even in the case of heterogenities and anisotropies.

The hybridizable discontinuous Galerkin method
Let T h be a regular triangulation ofΩ with elements K of diameter h K and let us define h := max K∈T h h K the mesh size. We denote by E h the set of all faces in the triangulation, and by E I and E Γ the set of interior and boundary faces, respectively. We use the standard notation for Sobolev spaces L 2 (Ω), H 1 0 (Ω), etc. We want to approximate the solution (u, p) of (1) with discrete functions where L 2 0 (Ω) := q ∈ L 2 (Ω) : Ω q = 0 and P r (K) is the set of polynomials in K with total degree less or equal than r.
To derive the hybridizable discontinuous Galerkin (HDG) formulation of (1) we introduce two new unknownsp andû. These new unknowns are usually called numerical fluxes and can be interpreted as single-valued approximations of p and u over E h , respectively.
For the numerical fluxes we use the following spaces defined only at the faces of the triangulation We write (·, ·) T h := K∈T h (·, ·) K and ·, · ∂T h := K∈T h ·, · ∂K where (·, ·) K and ·, · ∂K denotes the usual inner products in L 2 (K) and L 2 (∂K), respectively.
From now on we will assume that K| K is a continuous tensor and define

|36
Ingeniería y Ciencia thus by the Rayleigh-Ritz Theorem (a.k.a Min-Max theorem) we know that Finally, to obtain a weak formulation of problem (1) in the discrete spaces, we multiply it by convenient test functions and use Green's formulas to get: Here the expressions ∇ h and div h means the gradient and the divergence restricted to each element, respectively. The main goal of the hybridizable methods is to write the Problem P 1 G in terms of one of the numerical fluxes, let us sayp (since it has less number of degrees of freedom thanû). To that end we define the relation between the numerical fluxes given byû for some non-negative > 0, called stabilization parameter and defined over E h . Using (4) in Problem P 1 G and Green's identity, we get the following problem.
whereg is the extension by zero of the boundary condition g over E h . The last equation is introduced to impose the continuity of the normal trace ofû over the interior faces and the boundary condition.

Proof. It is enough to show that if
to the corresponding homogeneous problem, then it has to be the trivial (zero) solution. So let us consider f ≡ 0 andg ≡ 0. If we take v h = u h , q h = p h and µ =p in Problem P 2 G , we get after adding the resulting equations and simplifying, Finally, from the first equation of Problem P 2 G , using the Green's identity and the fact that Now that we have proved that Problem P 2 G has a unique solution, we can introduce the main advantage of the HDG method (over traditional DG methods) and reduce Problem P 2 G to a problem only over the skeleton. To this end we introduce the concept of local solver.

The local solvers
Notice that if the numerical fluxp is known, it is possible to find the solution In other words, if we considerp as a given data then (u h , p h ) ∈ V h × Q h satisfies the following local problem.
Problem P Loc can be written in terms of operators as and it is called local solver. It is easy to prove that the following result holds.
with A : are the solutions of the local solvers L(λ, 0) and L(0, f ), respectively. This theorem shows that the numerical fluxp can be found as the solution of a simpler problem (7) (since this problem is formulated only on the skeleton), and then the main unknowns (u h , p h ) can be recovered from the local solver Problem P Loc . This is the main feature of the HDG methods and shows how the degrees of freedom are reduced.

Error analysis
In this section we show the error analysis of the solution of Problem P 2 G based in projection operators on the approximation spaces Q h and V h . We follow the ideas of [25] given for the Laplace problem. Let The following lemma is a standard result and we refer to [26] for the proof.
Let us now define the projection errors given by where P W is the L 2 orthogonal projection over the space W h . These projection errors satisfy the problem given in the following lemma.
Proof. We just show the first equation, since the rest follow by similar arguments. First notice that the exact solution of (1) satisfies Then by the definition of the projectors in (8) and taking into account that divv h ∈ P r−1 (K) and the definition of P W , we get so if we subtract from (10) the first equation of Problem P 2 G we get Proof. This is a direct consequence of Lemma 3.2, taking v h = ε u h , q h = ε p h and µ = −εp h · n, adding the resulting equations.
Remark 3.1. Using Holder's inequality, inequality (2) and Corollary 3.1 we can find error estimates for ε u h in T h and ε p h − εp h on ∂T h ; however we still need to bound the error ε p h in T h . In order to accomplish this we use a classic duality argument. Let us define the auxiliary problem: Given Θ ∈ L 2 (Ω), find (σ, z) ∈ H(div; Ω) × H 1 m (Ω) such that This problem has a unique solution according with the Lax-Milgram theorem. Even more, if we assume that Ω is convex then we have that z ∈ H 2 (Ω) and Proof. The proof of this lemma is very similar to the proof of Lemma 4.1 in [25].
and taking supremum we get (13). For the inequality (14) we need to bound the right hand side of (13). First we take r = 0 in (9) (applied to the auxiliary problem) to obtain Finally, we complete the error analysis with the following theorem.

Ingeniería y Ciencia
We have proved that the HDG method is a super-convergent method for polynomial approximations of order r, since optimal convergence order is obtained for both p and u.
In the following section we present the high order polynomial basis that we will use in Section 5 to show the super convergence in the numerical tests.

Dubiner basis
To construct the spaces V h and Q h over triangles in R 2 or tetrahedrons in R 3 we use the Dubiner Basis. This basis is an orthogonal and complete set of functions that generates spaces of polynomials of order r > 0. The Dubiner basis were proposed in [23]. For triangular elements, the approximation space on the standard reference triangle is chosen as in [24].
Finally, for the space W h we just need an one-dimensional basis, so we choose the Legendre polynomials at the reference edge I := [−1, 1] (also to take advantage of their orthogonality property). These polynomials can be defined as φ j (x) := P 0,0 n (x), with cardinality d 1 − 1 and where j and n correspond to the bijection j = j(n, m) mentioned before.

Numerical examples
In this section we propose two numerical examples and test the behavior of the HDG method when one uses Dubiner basis for high order approximations. Our test cases are manufactured examples, here we construct the suitable source and boundary conditions for a given analytical solutions. In the first test case we consider an homogeneous porous media and in the second test case an heterogeneous domain.
In the test cases and for a fixed order of approximation r, the convergence orders respect to the refinement of the mesh are calculated as follow (2) and We also check the development of the HDG method with respect to the order of approximation. The convergence orders respect to the polynomial order r are calculated as follow we denote p r h the discontinuous approximate pressure p h when using a polynomial basis of order r.

Test case 1. Homogeneous case.
Consider an homogeneous domain Ω := [0, 1] 2 with constant permeability K = I (the identity matrix). The source term f , the flow at the boundary g and the flow velocity u are chosen such that the pressure is p(x, y) := sin (2πx) sin (2πy) .

Ingeniería y Ciencia
In Table 1 we show the data for the five triangulations used in the numerical test. 3.125E-02 4096 6080 Figure 1 shows the analytical and approximate solution of (1) and Problem P 2 G , respectively. The approximated solution in Figure 1 is computed using Dubiner basis of order five (r = 5) and over a mesh with 4096 elements. In Figure 1 we also show the comparison of the magnitude of the velocity computed analytically and using the same HDG method with r = 5.
Our aim is to show via numerical examples the convergence of the error proposed in the section 3 for the scalar and vectorial unknowns. Table 2 shows the history of convergence for the first five order approximations and using five different meshes. As expected, the history of convergence, in Table 2, shows that the solution p h converges to p and u h converges to u with order h r+1 .
We display in Figure 2 the convergence of the error for the higher order approximations. In Figure 2 the error is displayed in a log-log scale. The reference lines correspond to the expected orders of convergence, i.e. h 5 and h 6 respectively. Accordingly with the theory presented in Section 3, in this example we showed the super convergence of the HDG method for smooth solutions.
We finally test the behavior of the HDG method in terms of the approximation order r. As expected, the rates of convergence β L 2 and β H 1 goes to one when r increase and this result validate the HDG method as a p-method.
We take data so that p(x, y) := sin(xy) is the exact solution of (1) and in Figure 3 we show the permeability and the source function f . ing.cienc., vol. 16, no. 32, pp. 33-54, julio-diciembre. 2020.        Figure 4 shows the analytical and approximate solution of (1) and Problem P 2 G , respectively. Here we use the same mesh as in the first test case. In other words, the approximated solution in Figure 4 is computed using Dubiner basis with r = 5 and over a mesh with 4096 elements. In this test case we also use the meshes in Table 1. In the case of the heterogeneous media we also get the super-convergence of the HDG method. The history of convergence of the error for the second test case is shown in Table 4. As expected, the solution p h converges to p and u h converges to u with order h r+1 . This result proves the super-convergences of the HDG method that makes this method a suitable strategy to solve problems over more complex domains or with complex characteristics.
ing.cienc., vol. 16, no. 32, pp. 33-54, julio-diciembre. 2020.  In Figure 5 the convergence of the error for the highest order approximations is shown. We highlight that for the case r = 5 and the finest mesh the precision of the approximation reaches the precision of the numerical integration and the error computation is not accurate enough, however the firs three meshes allow us to obtain the predicted convergence rate.
Finally, the convergence of the error with respect to the order of the approximation is displayed in Table 5. With this we conclude that the HDG method presented could be also use as a p-method by fixing a mesh and increasing the order r, even in the case with heterogeneity in coefficients.

Conclusions
We have introduced the hibridizable discontinuous Galerking method (HDG) as a suitable idea for constructing high order approximations to the solution of Darcy's law. We proved the existence and uniqueness of the solutions of the global and local problems that defines the HDG method. Furthermore, we proved the error estimates for the numerical solutions by using projection operators.
In the numerical results we show the super-convergence of the HDG method. We considered two test cases with homogeneous and heterogeneous media and use Dubiner basis to construct high-order approximations. With this we showed that the HDG method is a suitable strategy to solve problems over heterogeneous domains, e.i. with non-constant coefficients.
We highlight that the application of such type of DG methods is vast. Further work goes in the direction of a-posteriori error estimates and real-life applications including computational efficiency tests. In terms of the applicability of this method we propose as future work the use of HDG methods with Dubiner basis on more complex domains, including random media and non-regular geometries. In other words, the extension of the numerical analysis in this paper can be extended to two-phase flow and real-life applications as fracture model simulations, CO 2 storage and environmental pollution.