A conservative hybrid method for Darcy flow

TL;DR

This paper introduces a hybrid mimetic spectral element method for Darcy flow that ensures mass conservation and inter-element continuity, resulting in sparse, mesh-invariant systems that are computationally efficient.

Contribution

It presents a novel hybrid spectral element formulation that maintains key physical and topological properties for Darcy flow, independent of mesh shape or size.

Findings

Produces extremely sparse algebraic systems

Efficiently assembled and solved per element

Invariant under mesh transformations

Abstract

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.

Tables

Table 1. Table 1: For 2D. Left: Number of total unknowns as a function of N 𝑁 N , for K = 3 × 3 𝐾 3 3 K=3\times 3 . Right: Number of unknowns as a function of K 𝐾 K , for N = 3 𝑁 3 N=3 .

$N$	Full system	$\lambda$ only	$\lambda$ / Full
$5$	$825$	$60$	0.07
$10$	$3000$	$120$	0.04
$15$	$6525$	$180$	0.03
$20$	$11400$	$240$	0.02
$25$	$17625$	$300$	0.02

Table 2. Table 2: For 3D. Left: Number of total unknowns as a function of N 𝑁 N , for K = 3 × 3 × 3 𝐾 3 3 3 K=3\times 3\times 3 . Right: Number of unknowns as a function of K 𝐾 K , for N = 3 𝑁 3 N=3 .

$N$	Full system	$\lambda$ only	$\lambda$ / Full
$5$	$16875$	$1350$	$0.08$
$10$	$121500$	$5400$	$0.04$
$15$	$394875$	$12150$	$0.03$
$20$	$918000$	$21600$	$0.02$
$25$	$1771875$	$33750$	$0.02$

Equations

$$\left\{\begin{array}[]{ll}\boldsymbol{u}+\mathbb{A}\ \nabla p&=0\\[4.73611pt] \nabla\cdot\boldsymbol{u}&=f\end{array}\right.\quad\text{in}\ \Omega\quad\text{and}\quad\left\{\begin{array}[]{lll}\Gamma&=\Gamma_{D}\cup\Gamma_{N}\\[4.73611pt] p&=\hat{p}&\text{on}\ \Gamma_{D}\\[4.73611pt] \boldsymbol{u}\cdot\boldsymbol{n}&=\hat{\boldsymbol{u}}_{\boldsymbol{n}}&\text{on}\ \Gamma_{N}\end{array}\right.\;,$$

$${\left(\mathbf{u},\mathbf{v}\right)}_{{\mathbb{A}}^{-1},\Omega }={\int }_{\Omega }\left(\mathbf{u},{\mathbb{A}}^{-1}\mathbf{v}\right)d\Omega ,$$

$$\begin{array}[]{ll}\mathcal{L}\left(\boldsymbol{u},p,\lambda;f\right)=&\frac{1}{2}\int_{\Omega_{K}}\boldsymbol{u}^{T}\mathbb{A}^{-1}\boldsymbol{u}\ \mathrm{d}\Omega_{K}+\int_{\Omega_{K}}p\left(\nabla\cdot\boldsymbol{u}-f\right)\mathrm{d}\Omega_{K}\\[4.73611pt] &+\int_{\partial\Omega_{K}\setminus\partial\Omega}\lambda\left(\boldsymbol{u}\cdot\boldsymbol{n}\right)\ \mathrm{d}\Gamma-\int_{\Gamma_{D}}\hat{p}\left(\boldsymbol{u}\cdot\boldsymbol{n}\right)\ \mathrm{d}\Gamma-\int_{\Gamma_{N}}\lambda\left(\hat{\boldsymbol{u}}\cdot\boldsymbol{n}\right)\ \mathrm{d}\Gamma\end{array}\;.$$

$$\left\{\begin{array}[]{lllll}\textcolor{black}{\left(\boldsymbol{v},\boldsymbol{u}\right)_{\mathbb{A}^{-1},\Omega_{K}}}+\left\langle\nabla\cdot\boldsymbol{v},p\right\rangle_{\Omega_{K}}+\left\langle\left(\boldsymbol{v}\cdot\boldsymbol{n}\right),\lambda\right\rangle_{\partial\Omega_{K}\setminus\partial\Omega}&=\left\langle\boldsymbol{v}\cdot\boldsymbol{n},\hat{p}\right\rangle_{\Gamma_{D}}&\quad\forall\ \boldsymbol{v}\in H(\mathrm{div};\Omega_{K})\\[4.73611pt] \left\langle q,\nabla\cdot\boldsymbol{u}\right\rangle_{\Omega_{K}}&=\left\langle q,f\right\rangle_{\Omega_{K}}&\quad\forall\ q\in L^{2}\left(\Omega_{K}\right)\\[4.73611pt] \left\langle\mu,\left(\boldsymbol{u}\cdot\boldsymbol{n}\right)\right\rangle_{\partial\Omega_{K}\setminus\partial\Omega}&=\left\langle\mu,\hat{\boldsymbol{u}}\cdot\boldsymbol{n}\right\rangle_{\Gamma_{N}}&\quad\forall\ \mu\in H^{\frac{1}{2}}\left(\partial\Omega_{K}\right)\end{array}\right.\;.$$

$$h_i\left(\xi \right)=\frac{\left({\xi }^2-1\right)L_N^{\prime }\left(\xi \right)}{N\left(N+1\right)L_N\left({\xi }_i\right)\left(\xi -{\xi }_i\right)},$$

$${\left(a^h,b^h\right)}_I={\vec{a}}^T{\mathbb{M}}^{(0)}\vec{b},\text{where}{\mathbb{M}}_{i,j}^{(0)}={\int }_{-1}^1{h}_i(\xi )h_j(\xi )\mathrm{d}\xi ,$$

$${⟨a^h,b^h⟩}_I={\tilde{\vec{a}}}^T\vec{b}:={\vec{a}}^T{\mathbb{M}}^{(0)}\vec{b}\Rightarrow \tilde{\vec{a}}={\mathbb{M}}^{(0)}\vec{a}.$$

$$e_j(\xi )=-\sum _{k=1}^{j-1}\frac{\text{d}{h}_k}{\text{d}\xi }(\xi ),\text{such that}{\int }_{{\xi }_{j-1}}^{{\xi }_j}{e}_i(\xi )={\delta }_{ij}.$$

$${\left(p^h,q^h\right)}_I={\vec{p}}^T{\mathbb{M}}^{(1)}\vec{q},\text{where}{\mathbb{M}}_{i,j}^{(1)}={\int }_{-1}^1{e}_i(\xi )e_j(\xi )\mathrm{d}\xi ,$$

$${⟨p^h,q^h⟩}_I={\tilde{\vec{p}}}^T\vec{q}:={\vec{p}}^T{\mathbb{M}}^{(1)}\vec{q}\Rightarrow \tilde{\vec{p}}={\mathbb{M}}^{(1)}\vec{p}.$$

$$\frac{d}{d\xi }{a}^h\left(\xi \right)=\frac{d}{d\xi }\sum _{i=0}^N{\mathsf{a}}_i{h}_i\left(\xi \right)=\sum _{i=1}^N\left({\mathsf{a}}_i-{\mathsf{a}}_{i-1}\right)e_i\left(\xi \right).$$

$${\mathbf{u}}^h\left(\xi ,\eta \right)=\sum _{i=0}^N\sum _{j=1}^N{{\mathsf{u}}_{\mathsf{x}}}_{i,j}{h}_i(\xi )e_j(\eta )\vec{\hat{ı}}+\sum _{i=1}^N\sum _{j=0}^N{{\mathsf{u}}_{\mathsf{y}}}_{i,j}{e}_i(\xi )h_j(\eta )\vec{\hat{ȷ}}.$$

$${\left({\mathbf{v}}^h,{\mathbf{u}}^h\right)}_{{\mathbb{A}}^{-1},{\Omega }_K}=\sum _K{\vec{v}}_{K_i}^T{\mathbb{M}}_{{\mathbb{A}}^{-1},K_i}^{(1)}{\vec{u}}_{K_i},$$

$$\mathbb{M}^{(1)}_{\mathbb{A}^{-1},{K_{i}}}=\int_{K_{i}}\left(\begin{array}[]{l}h_{i}(\textcolor{black}{x})\ e_{j}(\textcolor{black}{y})\\[4.73611pt] e_{i}(\textcolor{black}{x})\ h_{j}(\textcolor{black}{y})\end{array}\right)\times\mathbb{A}^{-1}\left(\textcolor{black}{x},\textcolor{black}{y}\right)\left(\begin{array}[]{l}h_{i}(\textcolor{black}{x})\ e_{j}(\textcolor{black}{y})\\[4.73611pt] e_{i}(\textcolor{black}{x})\ h_{j}(\textcolor{black}{y})\end{array}\right)\ \mathrm{d}{K_{i}}\;.$$

$$\begin{array}[]{ll}\nabla\cdot\boldsymbol{u}^{h}&=\frac{\partial}{\partial x}\sum_{i=0}^{N}\sum_{j=1}^{N}\mathsf{u_{x}}_{i,j}\ h_{i}(x)e_{j}(y)+\frac{\partial}{\partial y}\sum_{i=1}^{N}\sum_{j=0}^{N}\mathsf{u_{y}}_{i,j}\ e_{i}(x)h_{j}(y)\\[4.73611pt] &=\sum_{i,j=1}^{N}\left(\mathsf{u_{x}}_{i,j}-\mathsf{u_{x}}_{i-1,j}+\mathsf{u_{y}}_{i,j}-\mathsf{u_{y}}_{i,j-1}\right)e_{i}\left(x\right)e_{j}\left(y\right)\end{array}\;.$$

$${⟨q^h,\nabla \cdot {\mathbf{u}}^h⟩}_{{\Omega }_K}=\sum _K{\vec{q}}_{K_i}^T{\mathbb{E}}^{2,1}{\vec{u}}_{K_i},$$

$$\mathbb{E}^{2,1}=\left[\begin{array}[]{cccccccccccccccccccccccc}-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0\\ 0&0&0&-1&0&0&1&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&-1&0&0&1&0&0&0&0&-1&0&0&1&0&0&0&0&0&0\\ 0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0\\ 0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0\\ 0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&-1&0&0&1&0&0&0\\ 0&0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0&0\\ 0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&0&0&-1&0&0&1&0\\ 0&0&0&0&0&0&0&0&-1&0&0&1&0&0&0&0&0&0&0&0&-1&0&0&1\end{array}\right]\;.$$

$${⟨{\mu }^h,{\mathbf{u}}^h\cdot \mathbf{n}⟩}_{\partial {\Omega }_K\setminus \partial \Omega }=\sum _K{\vec{\mu }}_{K_i}^T\mathbb{N}{\vec{u}}_{K_i}={\vec{\mu }}^T{\mathbb{E}}_{\mathbb{N}}\vec{u},$$

$$\mathbb{E}_{\mathbb{N}}=\left[\begin{array}[]{lllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllllll}0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{1}&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&0&\scalebox{0.5}{-1}&0&0&0&0&0&0&0&0\end{array}\right]\;.$$

$$\left[\begin{array}[]{cc}\mathbb{A}&{\mathbb{E}_{\mathbb{N}}}^{T}\\[4.73611pt] \mathbb{E}_{\mathbb{N}}&0\end{array}\right]\left[\begin{array}[]{l}\vec{X}\\[4.73611pt] \vec{\lambda}\end{array}\right]=\left[\begin{array}[]{l}\vec{F}\\[4.73611pt] 0\end{array}\right]\;,$$

$$\mathbb{A}=\left[\begin{array}[]{llllllllllllllllllllll}\begin{array}[]{lll}\mathbb{M}^{(1)}_{\mathbb{A}^{-1},K_{1}}&{\mathbb{E}^{2,1}}^{T}\\[4.73611pt] \mathbb{E}^{2,1}&0\end{array}\\[4.73611pt] &&\begin{array}[]{lll}\mathbb{M}^{(1)}_{\mathbb{A}^{-1},K_{2}}&{\mathbb{E}^{2,1}}^{T}\\[4.73611pt] \mathbb{E}^{2,1}&0\end{array}\\[4.73611pt] &&&&\ddots&\ddots\\[4.73611pt] &&&&\ddots&\ddots\\[4.73611pt] &&&&&&\begin{array}[]{lll}\mathbb{M}^{(1)}_{\mathbb{A}^{-1},K}&{\mathbb{E}^{2,1}}^{T}\\[4.73611pt] \mathbb{E}^{2,1}&0\end{array}\end{array}\right]\;,$$

$$\vec{\lambda }={\left({\mathbb{E}}_{\mathbb{N}}{\mathbb{A}}^{-1}{{\mathbb{E}}_{\mathbb{N}}}^T\right)}^{-1}\cdot \left({\mathbb{E}}_{\mathbb{N}}{\mathbb{A}}^{-1}\vec{F}\right).$$

$$\left[\begin{array}[]{l}\vec{u}\\[4.73611pt] \vec{p}\end{array}\right]_{K_{i}}=\left[\begin{array}[]{ccc}\mathbb{M}^{(1)}_{\mathbb{A}^{-1}}&{\mathbb{E}^{2,1}}^{T}\\[4.73611pt] \mathbb{E}^{2,1}&0\end{array}\right]_{K_{i}}^{-1}\left[\begin{array}[]{l}\vec{\lambda}\\[4.73611pt] \vec{f}\end{array}\right]_{K_{i}}\;.$$

$$f_{ex}=\nabla \cdot (-\mathbb{A}\nabla p_{ex}),\text{where},$$

$$\begin{array}[]{ll}\mathbb{A}&=\frac{1}{x^{2}+y^{2}+\alpha}\left(\begin{array}[]{ll}10^{-3}x^{2}+y^{2}+\alpha&\quad\left(10^{-3}-1\right)xy\\[4.73611pt] \left(10^{-3}-1\right)xy&\quad x^{2}+10^{-3}y^{2}+\alpha\end{array}\right)\;;\quad\quad\alpha=0.1\\[4.73611pt] p_{ex}&=\mathrm{sin}\left(2\pi x\right)\mathrm{sin}\left(2\pi y\right)\end{array}\;,$$

Full text

1 Introduction

Hybrid formulations 2010Boffi ; 2015Cockburn ; 2018Zhang 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 2017Jain \textcolorblackwhich 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 \textcolorblackreduced 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 2008Herbin .

2 Darcy flow formulation

For $\Omega\in\mathbb{R}^{d}$, where $d$ is the dimension of the domain, the governing equations for Darcy flow, are given by,

[TABLE]

where, $\boldsymbol{u}$ is the velocity, $p$ is the pressure, $f$ the \textcolorblackprescribed source term, $\mathbb{A}$ is a $d\times d$ \textcolorblacksymmetric positive definite matrix, $\hat{p}$ and $\hat{\boldsymbol{u}}_{\boldsymbol{n}}$ are the \textcolorblackprescribed pressure and flux boundary conditions, respectively.

Notations

For $f,g\in L^{2}\left(\Omega\right)$, $\left(f,g\right)_{\Omega}$ denotes the usual $L^{2}$ - inner product.

For vector-valued function\textcolorblacks in $L^{2}$ we define the weighted inner product by,

[TABLE]

\textcolor

blackwhere $\left(\cdot\;,\cdot\right)$ denotes the pointwise inner product.

Duality pairing, denoted by $\left\langle\cdot,\cdot\right\rangle_{\Omega}$, is the outcome of a linear functional on $L^{2}\left(\Omega\right)$ acting on elements from $L^{2}\left(\Omega\right)$.

Let $\Omega_{K}$ be a \textcolorblackdisjoint partitioning of \textcolorblack$\Omega$ with total number of elements $K$, and $K_{i}$ is any element in $\Omega_{K}$, such that, $K_{i}\in\Omega_{K}$. We define the following broken Sobolev spaces 2016Carsten , \textcolorblack$H\left(\mathrm{div};\Omega_{K}\right)=\prod_{K}H\left(\mathrm{div};K_{i}\right)$, and $H^{1/2}\left(\partial\Omega_{K}\right)=\prod_{K}H^{1/2}\left(\partial K_{i}\right)$.

Weak formulation

The Lagrange functional for Darcy flow is defined as,

[TABLE]

The variational problem is then given by: For given \textcolorblack$f\in L^{2}\left(\Omega_{K}\right)$, find $\boldsymbol{u}\in H(\mathrm{div};\Omega_{K})$, \textcolorblack$p\in L^{2}\left(\Omega_{K}\right)$, \textcolorblack$\lambda\in H^{\frac{1}{2}}\left(\partial\Omega_{K}\right)$, such that,

[TABLE]

3 Basis functions

Primal and dual nodal degrees of freedom

Let $\xi_{j}$, $j=0,1,...,N$, be the $N+1$ Gauss-Lobatto-Legendre (GLL) points in $I\in\left[-1,1\right]$. The Lagrange polynomials $h_{i}(\xi)$ through $\xi_{j}$, of degree $N$, given by,

[TABLE]

form the 1D primal nodal polynomials which satisfy, $h_{i}(\xi_{j})=\delta_{ij}$.

Let $a^{h}$ and $b^{h}$ be two polynomials expanded in terms of $h_{i}\left(\xi\right)$. The $L^{2}$ - inner product is then given by,

[TABLE]

and, $\vec{a}=\left[\mathsf{a}_{0}\ \mathrm{a}_{1}\ \ldots\ \mathsf{a}_{N}\right]$ and $\vec{b}=\left[\mathsf{b}_{0}\ \mathsf{b}_{1}\ \ldots\ \mathsf{b}_{N}\right]$ are the nodal degrees of freedom. We define the algebraic dual degrees of freedom, $\widetilde{\vec{a}}$, such that the duality pairing is simply the vector dot product between primal and dual degrees of freedom,

[TABLE]

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 2011Gerritsma ,

[TABLE]

Let $p^{h}$ and $q^{h}$ be two polynomials expanded in edge basis functions. The inner product in $L^{2}$ space is given by,

[TABLE]

and, $\vec{p}=\left[\mathsf{p}_{1}\ \mathsf{p}_{2}\ \ldots\ \mathsf{p}_{N}\right]$ and $\vec{q}=\left[\mathsf{q}_{1}\ \mathsf{q}_{2}\ \ldots\ \mathsf{q}_{N}\right]$ are the edge degrees of freedom. As before, we define the dual degrees of freedom such that,

[TABLE]

\textcolor

blackA similar construction can be used for dual degrees of freedom in higher dimension. For construction of the dual degrees of freedom in $2D$ see 2017Jain and for $3D$ see 2018Zhang2 .

Differentiation of nodal polynomial representation

Let $a^{h}\left(\xi\right)$ be expanded in Lagrange polynomial, then

[TABLE]

Therefore, taking \textcolorblackthe derivative of \textcolorblacka polynomial involves two steps : First, \textcolorblacktake the difference of degrees of freedom; and second, change of basis from nodal to edge 2011Gerritsma .

4 Discrete inner product and duality pairing

For 2D domains, the higher dimensional primal basis are constructed \textcolorblackusing the tensor product of the 1D basis.

For the weak formulation in (2) we expand the velocity $\boldsymbol{u}^{h}$ in primal edge basis as,

[TABLE]

Weighted inner product

Using (1) and the expansions in (4), the weighted inner product is evaluated as,

[TABLE]

where, $\vec{u}_{K_{i}}$ are the degrees of freedom in element ${K_{i}}$, and

[TABLE]

For mapping of elements please refer to 2017Gerritsma .

Divergence of velocity

Divergence of velocity, $\nabla\cdot\boldsymbol{u}^{h}$, is evaluated using (3), but now for 2D,

[TABLE]

The pressure is expanded in the dual basis $\widetilde{e_{i}\left(\xi\right)e_{j}\left(\eta\right)}$. These basis are dual to the basis in which $\nabla\cdot\boldsymbol{u}^{h}$ is expanded in (5). Therefore the \textcolorblackweak constraint on divergence of velocity is a duality pairing evaluated as,

[TABLE]

where $\mathbb{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 \textcolorblackin $\Omega_{K}$. For an extensive discussion on the incidence matrix, see 2017Gerritsma . For an element of degree $N=3$,

[TABLE]

Connectivity matrix

\textcolor

blackThe connectivty matrix ensures continuity of the velocity across the elements. $\vec{\lambda}$ is the interface variable between the elements that acts as Lagrange multiplier that imposes the constraint given by,

[TABLE]

\textcolor

blackwhere $\mathbb{N}$ is the discrete trace operator. It is a sparse matrix that consists of $1$, $-1$ and [math] only. \textcolorblackFor construction of $\mathbb{N}$ please refer to 2018Gerritsma . \textcolorblack$\mathbb{E}_{\mathbb{N}}$ is the assembled $\mathbb{N}$ for all the elements. For discretization, $K=2\times 2$, $N=2$, $\mathbb{E}_{\mathbb{N}}$ is shown in (6). The matrix size of $\mathbb{E}_{\mathbb{N}}$ is $8\times 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.

[TABLE]

5 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,

[TABLE]

where, $\mathbb{A}$ is an invertible block diagonal matrix given by,

[TABLE]

$\mathbb{E}_{\mathbb{N}}$ is as given in (6), \vec{X}=\sum_{K}\left[\begin{array}[]{c}\vec{u}\\[4.73611pt] \vec{p}\end{array}\right]_{K_{i}}, and \vec{F}=\sum_{K}\left[\begin{array}[]{c}\vec{\hat{p}}\\[4.73611pt] \vec{f}\end{array}\right]_{K_{i}}, where $\vec{f}$ are the expansion coefficients of \textcolorblack$f^{h}\left(x,y\right)=\sum_{i,j}^{N}\mathsf{f}_{ij}\ e_{i}\left(x\right)e_{j}\left(y\right)$.

In (8), the mass matrix $\mathbb{M}^{(1)}_{\mathbb{A}^{-1},{K_{i}}}$ is the only dense matrix and also the only component that changes with each local element, ${K_{i}}$. $\mathbb{E}_{\mathbb{N}}$ is a sparse incidence matrix for the global system and $\mathbb{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 $\vec{\lambda}$, 2010Boffi ,

[TABLE]

\textcolor

blackTo evaluate $\vec{\lambda}$ in (9) we need $\mathbb{A}^{-1}$ that can be calculated efficiently by taking inverse of each block of $\mathbb{A}$ separately. This part can be easily parallelized. Once the $\vec{\lambda}$ is determined the solution in each element, ${K_{i}}$, can be evaluated by solving for,

[TABLE]

\textcolor

blackHere 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 \textcolorblackand easily parallelized.

\textcolor

blackThe 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 $\lambda$ 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 $\lambda$ 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 $\lambda$ system is certainly smaller than the full system, although the growth rate in size of $\lambda$ and full systems does not change significantly.

6 Results

In this section we present the results for a test problem from 2008Herbin by solving system (7). The domain of test problem is, $\Omega\in\left[0,1\right]^{2}$. The source term is defined as,

[TABLE]

[TABLE]

and Dirichlet boundary conditions are imposed along the entire boundary, $\Gamma_{\mathrm{D}}=\Gamma$ and $\Gamma_{\mathrm{N}}=\emptyset$.

We solve this problem on an orthogonal and a highly curved mesh, see Fig. 1.

The same problem was earlier addressed by authors in 2017Gerritsma , but for a method with continuous elements and primal basis only. For the configuration $K=3\times 3$, $N=6$, we compare the sparsity structure of the two approaches in Fig 2. On left we see the hybrid formulation, and on the right we see the continuous elements formulation \textcolorblack2017Gerritsma . The number of non zero entities are almost half in the hybrid formulation, $66384$, as 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 $\vec{\lambda}$ only system with continuous element system, \textcolorblackfor $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 \textcolorblackalmost $\mathcal{O}\left(10^{2}\right)$ higher. On the right we see the growth in condition number with increasing polynomial degree \textcolorblackfor $K=9\times 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 $\|\nabla\cdot u^{h}-f^{h}\|$. On the left side as a function of element size, $h=1/\sqrt{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 $\mathcal{O}\left(10^{-12}\right)$.

In Fig. 5, on the top two figures we show the error in the $H\left(\mathrm{div};\Omega\right)$ norm for the velocity; and at the bottom two figures we show the error in $L^{2}\left(\Omega\right)$ 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.