Efficient implementation of the hybridized Raviart-Thomas mixed method by converting flux subspaces into stabilizations

1.   Introduction

Let us begin by noting that the title of the paper contains both words: 'mixed methods' and 'stabilization'. At first glance, this might appear to be a technical error since mixed methods do not need stabilization since stability is directly ensured by the choice of their function spaces. However, this is not an error. In a recent paper, ^[5], one of the authors showed how a portion of the flux space of any hybridized mixed method can be recast as a stabilization of a equivalent hybridizable discontinuous Galerkin (HDG) methods. By equivalent, we mean that the original hybridized mixed method and the new HDG methods result in exactly the same solution. In this paper, we show that the implementation of the equivalent HDG method is faster than that of original hybridized mixed method. We carry out this for two equivalent HDG methods.

We apply this implementation procedure to a particular mixed method, namely, the Raviart-Thomas (RT) method ^[13] in symplexes for two reasons. The first is that its extension to any other mixed (or HDG) method defined on polytopal meshes is particularly straightforward^* The second reason is that, if this procedure does not work on a simple mixed method, it is very unlikely it would have a chance to work on more sophisticated mixed methods on polytopal meshes. In other words, we consider the material presented here as a necessary stepping stone towards the treatment of the cases of mixed and HDG methods defined in general polytopes.

^*How to define mixed methods for polyhedral elements has been shown in the series of papers on M-decompositions ^[6,7,9] and on new commuting diagrams ^[8]. See also the review ^[10] and the references therein.

The paper is organized as follows. In Section 2, we describe how to obtain an HDG method from a hybridized mixed method, we follow ^[5]. In Section 3, we describe in full detail the implementation of the hybridized RT method. We then do the same for two equivalent HDG methods. In Section 4, we display our numerical results. We end with some concluding remarks. We use the standard notation used for HDG methods, see, for example, ^[4,5].

2.   Background

2.1.   Subspace-to-stabilization.

For the sake of completeness, we begin by summarizing the subspace-to-stabilization result [5, Section 5]. Consider the Poisson problem:

with the boundary condition that $u = u_D$ on the boundary $\partial \Omega$. Here, $c$ is the coefficient, $f$ is the source term and $u_D$ is the Dirichlet boundary data. In this paper, we consider the case $c = Id$. The method can be easily generalized to arbitrary symmetric positive definite tensor fields $c$. The hybridized mixed method formulation for this problem is as follows: For each element $K$ of the mesh, find $\mathit{\boldsymbol{q}}_h$ and $u_h$ belonging to the local function spaces $\mathit{\boldsymbol{V}}(K)$ and $W(K)$, respectively, such that:

for all test functions $\mathit{\boldsymbol{v}}$ and $w$ belonging to the local function spaces $\mathit{\boldsymbol{V}}(K)$ and $W(K)$, respectively. The above equations are referred to as the 'local problem'. Here, $\hat{{u}}_h$ is an approximation to $u$ on the faces of the triangulation and is a data to the above problem. The additional equations for this face variable are:

where $F_I$ is each interior face of the mesh, $F_D$ is each Dirichlet boundary face of the mesh and $M_h(F)$ is the local function space on face $F$. The above equations are referred to as the 'global problem'.

Split $\mathit{\boldsymbol{V}}(K)$ into $\mathit{\boldsymbol{V}}_s(K) \oplus \mathit{\boldsymbol{V}}_a(K)$. Here, $\mathit{\boldsymbol{V}}_s(K)$ is the subspace of $\mathit{\boldsymbol{V}}(K)$ that will be converted into stabilization and $\mathit{\boldsymbol{V}}_a(K)$ is the subspace of $\mathit{\boldsymbol{V}}(K)$ that will be used to define the local problem of the equivalent HDG method. Their exact definition is deferred until later. This splitting converts the local problem (Eq (2.2)) into:

for all test functions $\mathit{\boldsymbol{v}}_a$, $\mathit{\boldsymbol{v}}_s$, and $w$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$, $\mathit{\boldsymbol{V}}_s(K)$, and $W(K)$, respectively. Requiring the functions in $\mathit{\boldsymbol{V}}_a(K)$ and $\mathit{\boldsymbol{V}}_s(K)$ to be orthogonal to each other in the $(\cdot, \cdot)_K$ inner-product, we obtain:

Integrating the second equation by parts and requiring $\mathit{\boldsymbol{V}}_s(K)$ to be any subspace of $\mathit{\boldsymbol{V}}(K)$ that is $L^2(K)$-orthogonal to $\mathit{\boldsymbol{\nabla}} W(K)$ yields

Note that in the above equation, $\mathit{\boldsymbol{q}}_s$ (which is the portion of $q$ in the subspace $\mathit{\boldsymbol{V}}_s(K)$) depends solely on the jump $u_h-\hat{{u}}_h$ on the faces of the element. Observe that the appearance of the term $u_h-\hat{{u}}_h$ has some similarity to the flux stabilization $\tau(u_h-\hat{{u}}_h)$ that is used in a HDG method (note that $(\tau(\cdot)$ here is a linear-mapping that satisfies certain requirements). This similarity is exploited to define a stabilization function $\tau(u_h-\hat{{u}}_h)$ for the HDG method that is equivalent to the above hybridized mixed method below.

Based on the form of Eq (2.5), let us define $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}$ to be the local lifting operator that maps a function $\mu$ in the space $L^2(\partial K)$ to the function $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(\mu)$ in the space $\mathit{\boldsymbol{V}}_s(K)$ as:

for all test functions $\mathit{\boldsymbol{v}}_s$ in the local space $\mathit{\boldsymbol{V}}_s(K)$. Then, $\mathit{\boldsymbol{q}}_s = \mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(u_h-\hat{{u}}_h)$, where by $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(u_h)$, we mean $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(u_h|_{\partial K})$. Moreover, the hybridized mixed method in Eq (2.4) can be manipulated to the following HDG method for the portion $\mathit{\boldsymbol{q}}_a$ of the flux approximation and the full scalar approximation $u_h$: Find $\mathit{\boldsymbol{q}}_a$ and $u_h$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively, such that:

where the stabilization function $\tau(\cdot)$ is defined using the lifting operator as: $\tau(u_h-\hat{{u}}_h)$ = $\mathit{\boldsymbol{n}}\cdot\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(u_h-\hat{{u}}_h)$. The global problem of the HDG method is same as the one in Eq (2.3).

The flux approximation $\mathit{\boldsymbol{q}}_h$ of the hybridized mixed method in Eq (2.2) is then obtained from the above HDG method as $\mathit{\boldsymbol{q}}_h = \mathit{\boldsymbol{q}}_a+\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}(u_h-\hat{{u}}_h)$. Thus, the effect of the local space $\mathit{\boldsymbol{V}}_s(K)$ is fully encapsulated in the defined stabilization function $\tau$ via the lifting operator $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}$. This is the crux of the 'spaces-to-stabilization' idea.

In summary, the conditions on the local spaces $\mathit{\boldsymbol{V}}_a(K)$ and $\mathit{\boldsymbol{V}}_s(K)$ are that:

The former and latter conditions above were used to obtain the Eqs (2.5) and (2.4), respectively.

2.2.   The equivalent hybridizable discontinuous Galerkin method after static condensation

Similar to other HDG methods, the degrees of freedom corresponding to $\mathit{\boldsymbol{q}}_a$ and $u_h$ can be statically condensed to yield a globally-coupled problem just for the degrees of freedom corresponding to the face variable $\hat{{u}}_h$ as follows. The local problem for the equivalent HDG method given in Eq (2.7) can be shown to be equal to the following local problem:

Find $\mathit{\boldsymbol{q}}_a$ and $u_h$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively, such that

for all test functions $\mathit{\boldsymbol{v}}_a$ and $w$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively. The influence of $\hat{{u}}_h$ and $f$ on $(\mathit{\boldsymbol{q}}_a, u_h)$ can be separated as $(\mathit{\boldsymbol{q}}_a, u_h)$ = $(\mathit{\boldsymbol{Q}}_{\hat{{u}}_h}, U_{\hat{{u}}_h})$ + $(\mathit{\boldsymbol{Q}}_f, U_f)$. Here, $(\mathit{\boldsymbol{Q}}_{\hat{{u}}_h}, U_{\hat{{u}}_h})$ is the solution to the following problem with $\mu = \hat{{u}}_h$:

Find $\mathit{\boldsymbol{Q}}_{\mu}$ and $U_{\mu}$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively, such that:

for all test functions $\mathit{\boldsymbol{v}}_a$ and $w$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively. $(\mathit{\boldsymbol{Q}}_f, U_f)$ is the solution to the problem:

Find $\mathit{\boldsymbol{Q}}_f$ and $U_f$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively, such that:

for all test functions $\mathit{\boldsymbol{v}}_a$ and $w$ in the local function spaces $\mathit{\boldsymbol{V}}_a(K)$ and $W(K)$, respectively.

Using the above decomposition, we can show that the equation for $\hat{{u}}_h$ given in Eq (2.3) is equal to the following problem:

Find $\hat{{u}}_h$ belonging to the global function space $M_h$ such that:

for each interior face $F_I$ and Dirichlet boundary face $F_D$ of the triangulation. Here, the global function space $M_h$ is the set of functions in $L^2(\mathcal{F})$, where $\mathcal{F}$ is the union of all faces $F$ in the mesh, and the restriction of $M_h$ on face $F$ is $M_h(F)$.

3.   Implementations of the hybridized Raviart-Thomas mixed method

We use this subspace-to-stabilization idea to come up with two new implementations of the hybridized RT mixed method. Each implementation stems from a different choice of the subspace $\mathit{\boldsymbol{V}}_s(K)$. We note that for the usual hybridized RT method, the local spaces $\mathit{\boldsymbol{V}}(K)$, $W(K)$ and $M_h(F)$ are:

Here, ${P}_k(K)$ denotes the space polynomials of degree $k$ defined on $K$, $\widetilde{{P}}_k(K)$ denotes the space of homogeneous polynomials of degree $k$ defined on $K$ and $d$ is the spatial dimension of the problem. Table 1 shows the three choices of $\mathit{\boldsymbol{V}}_s(K)$ and the name of the implementation that stems from each of these choices.

The implementation Usual-HRT is the usual implementation of the hybridized RT method, see ^[1]. This stems from the choice $\mathit{\boldsymbol{V}}_s(K) = \varnothing$ (the empty set). The implementations Stab-1-HRT and Stab-2-HRT are the two new implementations proposed in this paper. The new implementation Stab-1-HRT stems from the choice $\mathit{\boldsymbol{V}}_s(K) = \mathit{\boldsymbol{V}}_s^{(1)}(K)$, where $\mathit{\boldsymbol{V}}_s^{(1)}(K)$ is the largest subspace of $\mathit{\boldsymbol{V}}(K)$ that is $L^2-$orthogonal to $[P_k(K)]^d$. The other new implementation $\mathit{\boldsymbol{V}}_s^{(2)}(K)$ stems from choosing $\mathit{\boldsymbol{V}}_s(K) = \mathit{\boldsymbol{V}}_s^{(2)}(K)$, where $\mathit{\boldsymbol{V}}_s^{(2)}(K)$ is the space $\mathit{\boldsymbol{V}}_s^{(2)}(K)$ plus the vector-valued polynomials in $[P_k(K)]^d$ that are $L^2$-orthogonal to $[P_{k-1}(K)]^d$.

For each of these implementations, the local space $\mathit{\boldsymbol{V}}_a(K)$ is the $(\cdot, \cdot)_K$-orthogonal complement of $\mathit{\boldsymbol{V}}_s(K)$ within $\mathit{\boldsymbol{V}}(K)$. The space $\mathit{\boldsymbol{V}}_a(K)$ for the Usual-HRT, Stab-1-HRT, and Stab-2-HRT implementation is $[P_k(K)]^d+\mathit{\boldsymbol{x}}P_k(K)$, $[P_k(K)]^d$, and $[P_{k-1}(K)]^d$, respectively. The details of each of these implementations are given next.

3.1.   Usual-HRT (usual)

The details of the Usual-HRT implementation are given below.

3.1.1.   Basis

In the Usual-HRT implementation, the space $\mathit{\boldsymbol{V}}_a(K)$ equals $[P_k(K)]^d \oplus \mathit{\boldsymbol{x}}\widetilde{{P}}_k(K)$. We use the following basis functions for this space in each element $K$:

Here, $n = dm + m'$, where $m = C_{d}^{k+d-1}$, $m' = C_{d-1}^{k+d-1}$ and $C_a^b$ is the binomial coefficient $b!/(a!(b-a)!)$. The first $dm$ basis functions $\mathit{\boldsymbol{\varphi}}_1^{(K)}, \dots, \mathit{\boldsymbol{\varphi}}^{(K)}_{dm}$ correspond to the $[P_k(K)]^d$ portion of the space. The remaining $m'$ functions $\mathit{\boldsymbol{\varphi}}_{dm+1}^{(K)}, \dots, \mathit{\boldsymbol{\varphi}}_{n}^{(K)}$ correspond to the remaining portion of the space.

These basis functions are orthonormal to each other. They satisfy the orthonormality condition:

where $\delta_{ij}$ is the Kroenecker delta and $|K|$ is the measure (area in 2D and volume in 3D) of element $K$. The first $dm$ basis functions $\mathit{\boldsymbol{\varphi}}_1, \dots, \mathit{\boldsymbol{\varphi}}_{dm}$ are constructed using the (normalized) Dubiner polynomials ^[12] as:

Here, $\mathit{\boldsymbol{e}}_{j'}$ are the canonical basis functions of $\mathbb{R}^d$. $q_1, \dots, q_m$ are the orthonormal Dubiner polynomials in the element $K$ obtained by mapping the orthonormal polynomials in the reference simplex $\hat{{K}}$ to the element $K$ as:

where $\mathit{\boldsymbol{x}}^{(K)}(\hat{{\mathit{\boldsymbol{x}}}})$ is the affine mapping from the reference element $\hat{{K}}$ to the element $K$. These polynomials satisfy the following orthonormality relation:

The remaining $m'$ basis functions $\mathit{\boldsymbol{\varphi}}_{dm + 1}^{(K)}, \dots, \mathit{\boldsymbol{\varphi}}_n^{(K)}$ are constructed by multiplying the degree $k$ Dubiner polynomial basis functions with $\mathit{\boldsymbol{x}}$ and orthonormalizing them with the rest of the basis functions using the modified Gram-Schmidt kernel. This is described in Algorithm 1.

For the space $W(K)$ (which equals $P_k(K)$), we use the same orthonormal Dubiner polynomial basis functions:

For the space $M_h(F)$, we also use orthonormal Dubiner polynomial basis functions but defined on the faces of the mesh. They are:

where $F$ is a typical face of the mesh.

3.1.2.   Local problem

Using the above basis functions converts the local problem that depends on $\mu$ (Eq (2.10)) to the following matrix problem:

Here, $I_{n\times n}$ is the identity matrix, $0_{a\times b}$ is $a\times b$ matrix of zeros, $D^{(K)}_{m\times n}$ is the divergence matrix, $Q^{\mu, (K)}_{n\times (d+1)m'}$ is the degree of freedom matrix corresponding to the element-wise mapping $\mathit{\boldsymbol{Q}}_{\mu}$, $U^{\mu, (K)}_{m\times ((d+1)m')}$ is the degree of freedom matrix corresponding to the element-wise mapping $U_{\mu}$, and $-b^{Q^{\mu}, (K)}_{n\times ((d+1)m')}$ is the right-hand side matrix corresponding to the degrees of freedom $Q^{\mu, (K)}_{n\times (d+1)m'}$. The size of all matrices is given in the subscript. The entries of the divergence and right-hand matrices are:

respectively, where $F_j$ is the $j^{th}$ face of element $K$.

For efficient solution of the above local problem, we perform two optimizations. The first optimization pertains to computing the divergence matrix $D^{(K)}_{m\times n}$. This matrix has the following form:

Here, $D^{(K)^{(1)}}_{(m-m')\times dm}$ is the portion of the divergence-matrix from the first $dm$ basis functions, $\mathit{\boldsymbol{\varphi}}_1^{(K)}, \dots, \mathit{\boldsymbol{\varphi}}_{dm}^{(K)}$, which correspond to the $[P_k(K)]^d$ portion of the space. The remaining portion of the divergence matrix, $D^{(K)^{(2)}}_{m\times m'}$, is from the last $m'$ basis functions, $\mathit{\boldsymbol{\varphi}}_{dm+1}^{(K)}, \dots, \mathit{\boldsymbol{\varphi}}_{n}^{(K)}$. This form is exploited for its efficient computation as follows. The portion $D^{(K)^{(1)}}_{(m-m')\times dm}$ is first computed in along the coordinates of the reference element $\hat{{K}}$. We will denote this reference divergence matrix as $\hat{{D}}_{(m-m')\times dm}$ and is given by:

Then, to compute $D^{(K)^{(1)}}_{(m-m')\times dm}$ in each element $K$, we simply combine the columns of $\hat{{D}}_{(m-m')\times dm}$ using the entries of the Jacobian of $\mathit{\boldsymbol{x}}^{(K)}(\hat{{\mathit{\boldsymbol{x}}}})$. Hence, we do not differentiate the first $dm$ basis functions separately for each element of the mesh but just once in the reference element. Elementwise differentiation is performed only for the last $m'$ basis functions needed to compute the remaining of the divergence matrix, portion, $D^{(K)^{(2)}}_{m\times m'}$.

The second optimization pertains to the solution of the local matrix problem in Eq (3.4). Since the mass matrix is identity, the degrees of freedom matrix $Q^{\mu, (K)}_{n\times (d+1)m'}$ can be easily eliminated to obtain the following matrix problem just for the degrees of freedom $U^{\mu, (K)}_{m\times ((d+1)m')}$:

Here, $L^{(K)}_{m\times m}$ is the Laplacian matrix and is equal to:

To solve this problem, the matrix $L^{(K)}_{m\times m}$ is first factored using the (dense) Cholesky factorization method and the computed factor is used to obtain $U^{\mu, (K)}_{m\times ((d+1)m')}$.

Similarly, the above basis functions convert the local problem that depends on $f$ (Eq (2.11)) to the following matrix problem:

Here, $Q^{f, (K)}_{n\times 1}$ and $U^{f, (K)}_{m\times 1}$ are the degrees of freedom vector corresponding to the mapping $\mathit{\boldsymbol{Q}}_f$ and $U_f$, respectively.Here, $P^{(K)}f_{m\times 1}$ is the degree of freedom vector obtained by $L^2$-projection of $f$ onto $W$ using the above basis functions. Similar to the local problem that depended on $\mu$, the above matrix problem is also solved by eliminating the degrees of freedom corresponding to $Q^{f, (K)}_{n\times 1}$ and then using the above computed Cholesky factor to obtain $U^{f, (K)}_{m\times 1}$.

3.1.3.   Global problem

Using the above computed local problem solutions, the element matrix $A^{(K)}_{((d+1)m') \times ((d+1)m')}$ and element vector $b^{(K)}_{((d+1)m') \times 1}$ are computed in each element $K$ as:

respectively. These element matrices and vectors are assembled together and the degrees of freedom of $\hat{{u}}_h$ corresponding to the Dirichlet boundary faces are statically condensed to obtain the global matrix problem:

Here, $n_F$ is the total number of faces of the mesh minus the Dirichlet boundary faces. $A_{m'n_F \times m'n_{F}}$ is the global left-hand side matrix. $b_{m'n_{F} \times 1}$ is the global right-hand side vector. $\hat{{\mathrm{{u}}}}_{m'n_{F} \times 1}$ is the degree of freedom vector corresponding to $\hat{{u}}_h$. The above global matrix problem is solved using the sparse Cholesky factorization method.

3.2.   Stab-1-HRT (new)

The details of the new Stab-1-HRT implementation are given below. For general HDG methods, see ^[3,11].

3.2.1.   Basis

In this implementation, for the space $\mathit{\boldsymbol{V}}_a(K)$ in each element $K$, we use the following orthonormal basis functions:

The above functions $\mathit{\boldsymbol{\varphi}}_i^{(K)}$ are the same ones that were defined in Section 3.1.1. For the space $\mathit{\boldsymbol{V}}_s(K)$, we use:

These are the remaining $m'$ basis functions that were defined in Section 3.1.1. For the remaining spaces $W$ and $M_h$, we use the same basis functions as that used in Section 3.1.1.

3.2.2.   Local problem

In this implementation, we have to compute the stabilization mapping $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}$ (given in Eq (2.6)). Using the above basis functions for $\mathit{\boldsymbol{V}}_s(K)$, the matrix form of this mapping becomes:

Here, $b^{s, (K)}_{m'\times (d+1)m'}$ is the right-hand side of the stabilization mapping given by:

Using the above basis functions for $\mathit{\boldsymbol{V}}_a$ and the above matrix form of the stabilization mapping, the local problem that depends on $\mu$ (Eq (2.10)) becomes the following local matrix problem:

Here, $M^{s, (K)}_{m\times m}$ is the mass matrix that arises from the equivalent stabilization. It relates to the stabilization mapping matrix as:

where $B^{(K)}_{(d+1)m'\times m}$ is the linear mapping from the degrees of freedom of $u$ to that of $\mu$ in element $K$. $D^{(K)}_{m\times n}$ is the same divergence matrix that was defined in Section 3.1.2. The right-hand side matrix $b^{U^{\mu}, (K)}_{m\times ((d+1)m')}$ is given by:

Observe that the divergence matrix $D^{(K)}_{m\times dm}$ in Eq (3.10) has the following form:

where the portion $D^{(K)^{(1)}}_{(m-m')\times dm}$ is the same as that defined in Eq (3.6). Observe that the portion $D^{(K)^{(2)}}_{n\times m'}$ that was present in Eq (3.6) is not present above. Hence, we only compute the portion $D^{(K)^{(1)}}_{(m-m')\times dm}$ using the optimized implementation described in Section 3.1.2. We never compute $D^{(K)^{(2)}}_{n\times m'}$ that required differentiation of the last $m'$ basis functions separately in each element of the mesh. This yields in considerable reduction in the time taken to compute the individual local problems and is demonstrated by the numerical experiments in Section 4.

Similar to the Usual-HRT implementation, the degree of freedom matrix $Q^{\mu, (K)}_{dm\times (d+1)m'}$ is eliminated to obtain the following equation for $U^{\mu, (K)}_{m\times (d+1)m'}$:

Here, the Laplacian matrix $L^{(K), 1}_{m\times m}$ is given by:

The Laplacian is first factored using the Cholesky factorization method and the computed factors are used to compute $U^{\mu, (K)}_{m\times (d+1)m'}$.

Similarly, the chosen basis functions convert the local problem that depends on $f$ to the following matrix problem:

The above matrix problem is solved similar to the above local problem that depended only on $\mu$. The same Cholesky factors are used for its solution.

3.2.3.   Global problem

In this implementation, the element matrix and vector are computed as:

Then, similar to the Usual-HRT implementation, the above element matrices and element vectors are assembled to form the global problem:

This global problem is solved using the sparse Cholesky factorization.

3.3.   Stab-2-HRT (new)

The details of the new Stab-2-HRT implementation are given below.

3.3.1.   Basis

In this implementation, for the space $\mathit{\boldsymbol{V}}_a(K)$ in each element $K$, we use the following orthonormal basis functions:

where the above functions $\mathit{\boldsymbol{\varphi}}_i^{(K)}$ are the same functions defined in Section 3.1.1. For the space $\mathit{\boldsymbol{V}}_s(K)$, we use the remaining $(d+1)m'$ functions as the basis:

For the spaces $W$ and $M_h$, we use the same basis functions that were used in Section 3.1.1.

3.3.2.   Local problem

The local problem solution for this implementation is very similar to that of the Stab-1-HRT implementation. The matrix form of the stabilization mapping $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}$ is:

where

The local problem that depends on $\mu$ (Eq (2.10)) becomes the following matrix problem:

where the mass-matrix from the stabilization term is:

Then, similar to the Stab-1-HRT implementation, the above local matrix problem is solved using the Cholesky factorization methodology after eliminating the matrix $Q^{\mu, (K)}_{d(m-m')\times (d+1)m'}$.

The local problem that depends on $f$ becomes:

The above matrix problem is also solved using the same Cholesky factors that were computed while solving the above local problem that depended on $\mu$ alone.

3.3.3.   Global problem

In this implementation, the element matrix and vector are computed as:

Then, similar to the Stab-1-HRT implementation, the above element matrices and element vectors are assembled to form the global problem

and this global problem is solved using the sparse Cholesky factorization.

4.   Numerical results

We present results comparing the three implementations. We consider the Poisson problem with $f = 8\pi^2 \sin(2\pi x_1)\sin(2\pi x_2)$ in the domain $(0, 1)^2$. The domain is first split into 16 uniform quadrilateral elements along each direction. Each quadrilateral element is further split into two triangular elements. This leads to a uniform mesh of 512 triangular elements. Polynomial degrees 1 to 20 are considered. Zero-Dirichlet boundary condition is imposed on all the four sides of the domain. All these implementations were first validated to make sure they yield identical solutions.

All numerical experiments were performed in MATLAB. We use a workstation with Intel(R) Core(TM) i7-8700 processor. The processor has six cores and hyperthreading. We also note that MATLAB uses the multi-threaded MKL BLAS backend for certain matrix and vector manipulations. Hence, there is some inherent parallelism in our implementation. For the global problem solution, we use the sparse direct solver available in MATLAB. MATLAB uses the sparse multi-threaded Cholesky solver CHOLMOD ^[2] for our symmteric positive definite problem when performing $\texttt{x = A \\ b}$.

Tables 2 and 3 show the time consumed by the one-time operations (that are performed just once in the reference element), the local problem solution in all the elements, global problem solution and the total solution time. The breakdown of the time taken by the different components of the local problem solution are shown in Tables 4 and 5. The percentage benefit of the new implementations compared to the Usual-HRT implementation is shown in Table 6.

The new implementations Stab-1-HRT and Stab-2-HRT are faster than the Usual-HRT implementation for all the polynomial degrees. From Table 6, we observe that they are 10-20% faster depending on the polynomial degree. Stab-1-HRT is slightly faster than Stab-2-HRT for nearly all polynomial degrees (except degree six).

The local problem solutions consume the majority of the total solution time. These local problem solution solutions are faster for the new Stab-1-HRT and Stab-2-HRT implementations compared to the usual-HRT implementations. From Tables 4 and 5 (which show the breakdown of the local problem solutions), we observe that this performance benefit is essentially because in the new Stab-1-HRT and Stab-2-HRT implementations, we need not compute the derivative of the additional RT basis functions, i.e., $\mathit{\boldsymbol{\varphi}}_{dm+1}$, $\dots$, $\mathit{\boldsymbol{\varphi}}_n$, and their contribution to the divergence matrix of the local problems. However, in the usual-HRT implementation, these derivatives and their contribution to the divergence matrix must be computed. In the Stab-1-HRT and Stab-2-HRT implementations, there is an overhead associated with the computation of the stabilization mappings. However, the benefits from not computing the derivative of the additional RT basis functions and their contribution to the divergence matrix significantly outweights this overhead. Hence, the new implementations of the hybridized RT method yield significant (10–20%) performance benefit compared to the usual implementation.

5.   Conclusions and ongoing work

As pointed out in ^[5], although the choice of the local function space $\mathit{\boldsymbol{V}}_a(K)$ is not unique, as we have also seen here, the smallest of these local spaces, $\nabla W(K) = \nabla {P}_k(K)$, is actually unique. It remains to be seen it we still retain an advantage over the usual implementation of the hybridized RT method for this choice. This constitutes the subject of ongoing work.

The present paper is currently being extended along the following directions:

(1) As pointed out in the Introduction, the choice of the RT method as the mixed method is by no means restrictive, as the approach proposed here can also be applied to any other mixed method for polyhedral meshes, for e.g., those defined using M-decompositions in ^[6,7,9] and using new commuting diagrams ^[8] (see also the review ^[10]) and the references therein). The space $\mathit{\boldsymbol{V}}(K)$ in these methods also have the form $[P_k(K)]^d\oplus \mathit{\boldsymbol{V}}_{\text{fill}}$ and the $W(K)$ is still equal to $P_k(K)$. Similar to the implementation presented here, the effect of $\mathit{\boldsymbol{V}}_{\text{fill}}$ (and also a portion of $[P_k(K)]^d$) can be encapsulated within a mapping $\mathit{\boldsymbol{L}}_{\mathit{\boldsymbol{V}}_s}$ defined in a similar way.

In our new implementations, we have not defined the basis functions via the Piola Transform. On the other hand, the reference unit simplex is only used to define the orthonormal polynomial basis functions for $[P_k(K)]^d$ and to compute the divergence matrix via chain rule. In the case of polytopal elements, we note that the polytopal mixed methods in ^[6,7,8,9,10] do not make use of a reference element. Hence, it would be a different baseline to compare our implementations to.

Furthermore, the application of this approach to other general HDG methods can be carried out very easily, as we are going to show elsewhere.

(2) In $d = 3$ dimensions, we expect similar or better speedups. For e.g., for simplexes and Stab-2-HRT, the reduction in the computational effort for the local problem is approximately proportional to $dim(\mathit{\boldsymbol{V}}_s^{(2)}(K))$ which is $(d+1)\, C_{d-1}^{k+d-1}$. The factor

For large $k$, this factor is around $4/3$ times larger for $d = 3$ compared to $d = 2$. Hence, we expect similar or bigger speedups for three-dimensional problems. Extension to $d = 3$ is part of our ongoing work.

(3) Finally, note that, although we have exploited the fact that the tensor-valued function $c$ is the identity for our model second-order elliptic problem, it is easy to extend what has been done to a general elliptic problem. This is also part of our ongoing work.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

Bernardo Cockburn's research was supported in part by the Advanced Computational Center for Entry Systems Simulation (ACCESS) through NASA grant 80NSSC21K1117.

Conflict of interest

The authors declare no conflicts of interest.