deal.II Implementation of a Weak Galerkin Finite Element Solver for Darcy Flow

Abstract

This paper presents a weak Galerkin (WG) finite element solver for Darcy flow and its implementation on the deal.II platform. The solver works for quadrilateral and hexahedral meshes in a unified way. It approximates pressure by Q-type degree $k({\ge }0)$ polynomials separately defined in element interiors and on edges/faces. Numerical velocity is obtained in the unmapped Raviart-Thomas space $RT_{[k]}$ via postprocessing based on the novel concepts of discrete weak gradients. The solver is locally mass-conservative and produces continuous normal fluxes. The implementation in deal.II allows polynomial degrees up to 5. Numerical experiments show that our new WG solver performs better than the classical mixed finite element methods.

Harper, Liu, and Wang were partially supported by US National Science Foundation grant DMS-1819252. We thank Dr. Wolfgang Bangerth for the computing resources.

1 Introduction

The Darcy equation, although simple, plays an important role for modeling flow in porous media. The equation usually takes the following form

(1)

where $\varOmega$ is a 2-dim or 3-dim bounded domain, p is the unknown pressure, $\mathbf {K}$ is a conductivity matrix that is uniformly symmetric positive definite (SPD), c is a known function, f is a known source term, $p_D$ is a Dirichlet boundary condition, $u_N$ is a Neumann boundary condition, and $\mathbf {n}$ is the outward unit normal vector on $\partial \varOmega$, which has a nonoverlapping decomposition $\varGamma ^D \cup \varGamma ^N$.

The elliptic boundary value problem (1) can be solved by many types of finite element methods. But in the context of Darcy flow, local mass conservation and normal flux continuity are two most important properties to be respected by finite element solvers.

• The continuous Galerkin (CG) methods [5] use the least degrees of freedom but do not possess these two properties and hence cannot be used directly. Several post-processing procedures have been developed [7, 8].

• Discontinuous Galerkin (DG) methods are locally conservative by design and gain normal flux continuity after post-processing [4].

• The enhanced Galerkin (EG) methods [19] possess both properties but need to handle some minor issues in implementation.

• The mixed finite element methods (MFEMs) [2, 6] have both properties by design but result in indefinite discrete linear systems, for which hybridization needs to be employed to convert them into definite linear systems.

• The weak Galerkin (WG) methods [11, 13, 15,16,17, 20] have both properties and result in SPD linear systems that are easier to solve.

In this paper, we investigate efficient implementation of WG Darcy solvers in deal.II, a popular finite element package [3], with the intention to make WG finite element methods practically useful for large-scale scientific computation.

2 A WG Finite Element Solver for Darcy Flow

WG solvers can be developed for Darcy flow on simplicial, quadrilateral or hexahedral, and more general polygonal or polyhedral meshes. These finite element schemes may or may not contain a stabilization term, depending on choices of the approximating polynomials for pressure in element interiors and on edges/faces. Through integration by parts, these polynomial basis functions are used for computing discrete weak gradients, which are used to approximate the classical gradient in the variational form for the Darcy equation. Discrete weak gradients can be established in a general vector polynomial space [18] or a specific one like the Raviart-Thomas space [11, 17] that has desired approximation properties.

This paper focuses on quadrilateral and hexahedral meshes, in which faces are or very close to being flat. We use $Q_k (k \ge 0)$-type polynomials in element interiors and on edges/faces for approximating the primal variable pressure. Their discrete weak gradients are established in local unmapped Raviart-Thomas $RT_{[k]} (k \ge 0)$ spaces, for which we do not use the Piola transformation. We use the same form of polynomials as that for rectangles and bricks in the classical MFEMs [6].

To illustrate these new ideas, we consider a quadrilateral E centered at $(x_c,y_c)$. We define the local unmapped Raviart-Thomas space $RT_{[0]}(E)$ as

$$\begin{aligned} RT_{[0]}(E) = \text{ Span }(\mathbf {w}_1, \mathbf {w}_2, \mathbf {w}_3, \mathbf {w}_4), \end{aligned}$$

(2)

where

$$\begin{aligned} \mathbf {w}_1 = \left[ \begin{array}{c} 1 \\ 0 \end{array} \right] , \quad \mathbf {w}_2 = \left[ \begin{array}{c} 0 \\ 1 \end{array} \right] , \quad \mathbf {w}_3 = \left[ \begin{array}{c} X \\ 0 \end{array} \right] , \quad \mathbf {w}_4 = \left[ \begin{array}{c} 0 \\ Y \end{array} \right] , \end{aligned}$$

(3)

and $X=x-x_c$, $Y=y-y_c$ are the normalized coordinates.

Now we introduce a new concept of 5 discrete weak functions $\phi _i (0 \le i \le 4)$.

• $\phi _0$ is for element interior: It takes value 1 in the interior $E^\circ$ but 0 on the boundary $E^\partial$ (all 4 edges);

• $\phi _1, \phi _3, \phi _3, \phi _4$ are for the four sides respectively: $\phi _i (1 \le i \le 4)$ takes value 1 on the i-th edge but 0 on all other three edges and in the interior.

Any such function $\phi$ has two independent parts: $\phi ^\circ$ is defined in $E^\circ$, whereas $\phi ^\partial$ is defined on $E^\partial$, together written as $\phi = \{ \phi ^\circ , \phi ^\partial \}$. Its discrete weak gradient $\nabla _w \phi$ can be specified in $RT_{[0]}(E)$ via integration by parts [20]:

$$\begin{aligned} \displaystyle \int _{E} (\nabla _w \phi ) \cdot \mathbf {w} = \int _{E^\partial } \phi ^\partial (\mathbf {w} \cdot \mathbf {n}) - \int _{E^\circ } \phi ^\circ (\nabla \cdot \mathbf {w}), \quad \quad \forall \mathbf {w} \in RT_{[0]}(E). \end{aligned}$$

(4)

This attributes to solving a size-4 SPD linear system. Note that

• (i) For a quadrilateral, $\phi ^\circ$ or $\phi ^\partial$ each can also be a degree $k \ge 1$ polynomial and the discrete weak gradient $\nabla _w \phi$ is then established in the local unmapped Raviart-Thomas space $RT_{[k]} (k \ge 1)$.

• (ii) For a hexahedron with nonflat faces, we can use the averaged normal vectors in (4). The Jacobian determinant is still used in computation of the integrals.

For a rectangle $E = [x_1,x_2] \times [y_1,y_2]$ ($\varDelta x = x_2 - x_1, \varDelta y = y_2 - y_1$), we have

$$\begin{aligned} \left\{ \begin{array}{rrrrr} \displaystyle \nabla _w \phi _0 = &{} 0 \mathbf {w}_1 \;\; + &{} 0 \mathbf {w}_2 \;\; + &{} \frac{-12}{(\varDelta x)^2} \mathbf {w}_3 \;\; + &{} \frac{-12}{(\varDelta y)^2} \mathbf {w}_4, \\ \displaystyle \nabla _w \phi _1 = &{} \frac{-1}{\varDelta x} \mathbf {w}_1 \;\; + &{} 0 \mathbf {w}_2 \;\; + &{} \frac{6}{(\varDelta x)^2} \mathbf {w}_3 \;\; + &{} 0 \mathbf {w}_4, \\ \displaystyle \nabla _w \phi _2 = &{} \frac{1}{\varDelta x} \mathbf {w}_1 \;\; + &{} 0 \mathbf {w}_2 \;\; + &{} \frac{6}{(\varDelta x)^2} \mathbf {w}_3 \;\; + &{} 0 \mathbf {w}_4, \\ \displaystyle \nabla _w \phi _3 = &{} 0 \mathbf {w}_1 \;\; + &{} \frac{-1}{\varDelta y} \mathbf {w}_2 \;\; + &{} 0 \mathbf {w}_3 \;\; + &{} \frac{6}{(\varDelta y)^2} \mathbf {w}_4, \\ \displaystyle \nabla _w \phi _4 = &{} 0 \mathbf {w}_1 \;\; + &{} \frac{1}{\varDelta y} \mathbf {w}_2 \;\; + &{} 0 \mathbf {w}_3 \;\; + &{} \frac{6}{(\varDelta y)^2} \mathbf {w}_4. \end{array} \right. \end{aligned}$$

(5)

Let $\mathcal {E}_h$ be a shape-regular quadrilateral mesh. Let $\varGamma _h^D$ be the set of all edges on the Dirichlet boundary $\varGamma ^D$ and $\varGamma _h^N$ be the set of all edges on the Neumann boundary $\varGamma ^N$. Let $S_h$ be the space of discrete shape functions on $\mathcal {E}_h$ that are degree k polynomials in element interiors and also degree k polynomials on edges. Let $S_h^0$ be the subspace of functions in $S_h$ that vanish on $\varGamma ^D_h$. For (1), we seek $p_h = \{p_h^\circ ,p_h^\partial \} \in S_h$ such that $p_h^\partial |_{\varGamma ^D_h} = Q_h^\partial (p_D)$ (the $L^2$-projection of Dirichlet boundary data into the space of degree k polynomials on $\varGamma ^D_h$) and

$$\begin{aligned} \mathcal {A}_h(p_h,q) = \mathcal {F}(q), \qquad \forall q=\{q^\circ ,q^\partial \} \in S_h^0, \end{aligned}$$

(6)

where

$$\begin{aligned} \displaystyle \mathcal {A}_h(p_h,q) = \sum _{E \in \mathcal {E}_h} \int _E \mathbf {K} \nabla _w p_h \cdot \nabla _w q + \sum _{E \in \mathcal {E}_h} \int _E c \, p \, q, \end{aligned}$$

(7)

$$\begin{aligned} \displaystyle \mathcal {F}(q) = \sum _{E \in \mathcal {E}_h} \int _E f q^\circ - \sum _{\gamma \in \varGamma _h^N} \int _\gamma u_N q^\partial . \end{aligned}$$

(8)

This results in a symmetric positive-definite discrete linear system [17].

Note $\nabla _w p_h$ is in the local Raviart-Thomas space, but $-\mathbf {K} \nabla _w p_h$ may not be. A local $L^2$-projection $\mathbf {Q}_h$ is needed [11, 13, 17] to get it back into the RT space:

$$\begin{aligned} \mathbf {u}_h = \mathbf {Q}_h (-\mathbf {K} \nabla _w p_h). \end{aligned}$$

(9)

This is the numerical Darcy velocity for subsequent applications, e.g., transport simulations. Clearly, this process is readily parallelizable for large-scale computation. This numerical velocity is locally mass-conservative and the corresponding normal flux is continuous across edges or faces, as proved in [11, 17].

As shown in [17], this Darcy solver is easy to be implemented and results in a symmetric positive-definite system that can be easily solved by a conjugate-gradient type linear solver. The WG methodology has connections to but is indeed different than the classical mixed finite element methods, especially the hybridized MFEMs [13, 14].

3 deal.II Implementation of WG Solver for Darcy Flow

deal.II is a popular C++ finite element package [3]. It uses quadrilateral and hexahedral meshes instead of simplicial meshes. The former may involve less degrees of freedom than the latter. The resulting linear systems may have smaller sizes, although the setup time for these linear systems may be longer. The setup time is spent on bilinear/trilinear mappings from the reference square/cube to general quadrilaterals/hexahedra and computation of various integrals.

3.1 Quadrilateral and Hexahedral Meshes

deal.II handles meshes by the GridGenerator class. All mesh information, such as the number of active cells, degrees of freedom, are stored in this class. For any integer $k \ge 0$, our WG$(Q_k,Q_k;RT_{[k]})$ solver is locally mass-conservative and produces continuous normal fluxes regardless of the quality of quadrilateral and hexahedral meshes. In order to obtain the desired order k convergence rate in pressure, velocity, and normal fluxes, we require meshes to be asymptotically parallelogram or parallelopiped [11, 17].

3.2 Finite Element Spaces

The WG$(Q_k, Q_k; RT_{[k]})$ solver involves three finite element spaces. The first two spaces are for the pressure unknowns, the third one (RT space) is used for discrete weak gradients and numerical velocity. In deal.II implementation, the first two are combined as

The third one (with dim being 2 or 3) is

Raviart-Thomas Spaces for Discrete Weak Gradients and Velocity. WG allows use of unmapped RT spaces on quadrilaterals and hexahedra [11, 17]. These spaces use the same polynomials for shape functions as those in the classical RT spaces for 2-dim or 3-dim rectangles [6]. They are respectively,

$$\begin{aligned} RT_{[k]}(E) = Q_{k+1,k} \times Q_{k,k+1}, \end{aligned}$$

(10)

$$\begin{aligned} RT_{[k]}(E) = Q_{k+1,k,k} \times Q_{k,k+1,k}\times Q_{k,k,k+1}. \end{aligned}$$

(11)

In deal.II, we use degree for k in Eqs. (10) or (11) and have

Two Separate Polynomial Spaces for Pressure. Note that for the WG$(Q_k, Q_k; RT_{[k]})$ finite element solver for Darcy flow, the pressure is approximated separately in element interiors by $Q_k$-type polynomials and on edges/faces by $Q_k$-type polynomials also. Note that the 2nd group of $Q_k$-type polynomials are defined locally on each edge/face. For the deal.II implementation, we have

where degree is k, that is, the degree of the polynomials, “1” means these two groups of pressure unknowns are just scalars. Note that

• FE_DGQ is a finite element class in deal.II that has no continuity across faces or vertices, i.e., every shape function lives exactly in one cell. So we use it to approximate the pressure in element interiors.

• FE_FaceQ is a finite element class that is defined only on edges/faces.

However, these two different finite element spaces are combined into one finite element system, we split these shape functions as

Here “0” corresponds to the 1st finite element class FE_DGQ for the interior pressure; “1” corresponds to the 2nd finite element class FE_FaceQ for the face pressure. Later on, we will just use fe_values[interior].value and fe_values[face].value for assembling the element-level matrices.

3.3 Gaussian Quadratures

Finite element computation involves various types of integrals, which are discretized via quadratures, e.g., Gaussian quadratures. For example, we consider

$$\begin{aligned} \int _E f \approx \sum _{k=1}^{K} w_k f(x_k,y_k) J_k, \end{aligned}$$

(12)

where K is the number of quadrature points, $(x_k,y_k)$ is the k-th quadrature point, $J_k$ is the corresponding Jacobian determinant, and $w_k$ is the weight. In deal.II, this is handled by the Quadrature class. In particular, the Jacobian determinant value and weight for each quadrature point are bundled together as

where $q_k$ is the k-th quadrature point.

3.4 Linear Solvers

deal.II provides a variety of linear solvers that are inherited from PETSc. The global discrete linear systems obtained from the weak Galerkin finite element discretization of the Darcy equation are symmetric positive-definite. Thus we can choose a conjugate-gradient type linear solver for them.

3.5 Graphics Output

In our WG$(Q_k, Q_k; RT_{[k]})$ solver for Darcy flow, the scalar pressures are defined separately in element interiors and on edges/faces of a mesh. These values are output separately in deal.II. The interior pressures are handled by DataOut, whereas the face pressures are handled by DataOutFace. Specifically,

are used to subdivide each cell into smaller patches, which provide better visualization if we use higher degree polynomials. The post-processed data are saved as vtk files for later visualization in VisIt.

4 Code Excerpts with Comments

This section provides some code excerpts with comments. More details can be found in deal.II tutorial Step-61 (subject to minor changes) [1].

4.1 Construction of Finite Element Spaces

Note that FE_RaviartThomas is a Raviart-Thomas space for vector-valued functions, FESystem defines finite element spaces in the interiors and on edges/faces. Shown below is the code for the lowest order WG finite elements.

4.2 System Setup

The following piece distributes degrees of freedom for finite element spaces.

The following piece sets up matrices and vectors in the system.

4.3 System Assembly

The following piece uses extractors to extract components of finite element shape functions.

The following pieces calculates the Gram matrix for the RT space.

The following piece handles construction of WG local matrices.

The following piece calculates the local right-hand side.

The following piece distributes entries of local matrices into the system matrix and also incorporates the local right-hand side into the system right-hand side.

5 Numerical Experiments

This section presents three numerical examples (Eq. (1) with $c=0$) to demonstrate accuracy and robustness of our novel WG solver for Darcy flow.

Table 1. Ex.1: Convergence rates of WG$(Q_k,Q_k;RT_{[k]})$ solver on rectangular meshes

Example 1

(A smooth example for convergence rates). Here we have domain $\varOmega = (0,1)^2$, conductivity $\mathbf {K} = \mathbf {I}_2$, and a known solution for the pressure:

$$\begin{aligned} p(x,y) = \sin (\pi x)\sin (\pi y). \end{aligned}$$

A homogeneous Dirichlet boundary condition is posed on the entire boundary.

Fig. 1.

Ex.1: Numerical pressure by WG$(Q_1,Q_1;RT_{[1]})$ solver on rectangular meshes

The WG$(Q_k, Q_k; RT_{[k]})$ solver is tested on Example 1 for $k=0,1,2$ on a sequence of uniform rectangular meshes. As shown in Table 1, the solver exhibits order k convergence rates for the $L^2$-norms of the errors in the interior pressure, velocity, and normal flux. Shown in Fig. 1 are the profiles of the numerical pressure obtained from applying the WG$(Q_1, Q_1; RT_{[1]})$ solver. In the right panel, the edge pressures are plotted as grey line segments. The graphical results in both panels demonstrate nice monotonicity in the numerical pressure produced by our WG solver.

Fig. 2.

Example 2 (Heterogeneous permeability): Numerical pressure and velocity

Table 2. Example 2: Comparison between WG and MFEM solvers

Example 2

(Heterogeneous permeability). The permeability profile is adopted from [9]. We consider a simple Darcy flow problem on the unit square. Dirichlet boundary conditions are posed on the left and right sides: $p=1$ for $x=0$; and $p=0$ for $x=1$. The other two sides have a homogeneous Neumann (no-flow) boundary condition. The problem was also tested using Matlab in [17].

Shown in Fig. 2 right panel are the numerical pressure and velocity profiles obtained from apply our WG$(Q_0,Q_0;RT_{0})$ solver on a uniform $40 \times 40$ rectangular mesh. Clearly, the elementwise numerical pressure stays between 0 and 1, the pressure profile demonstrates monotonicity from left to right, and the velocity profile reveals the low-permeability regions and channels for fast flow.

Example 2 was also solved by a mixed finite element solver built in deal.II that is based on Schur complement (See deal.II tutorial Step-20). We compare the lowest order WG solver ($k=0$) with the lowest order MFEM solver on a sequence of rectangular meshes on a Toshiba laptop. The tolerance for linear solvers is $10^{-9}$. Table 2 shows that the WG solver produces very close results with significantly less runtime.

Fig. 3.

Example 3 (SPE10 Model 2): Permeability profiles on $log_{10}$ scales

Example 3

(Permeability profile in SPE10 Model 2). SPE10 was developed as a benchmark for upscaling methods, but the 2nd dataset is becoming a popular testcase for comparing different numerical methods. The dataset is a 3-dim geo-statistical realization from the Jurassic Upper Brent formations [12]. The model has geometric dimensions 1200 (ft) $\times$ 2200 (ft) $\times$ 170 (ft). The dataset is provided on a 60 $\times$ 220 $\times$ 85 Cartesian grid, in which each block has a size 20 (ft) $\times$ 10 (ft) $\times$ 2 (ft). The top 70 ft (35 layers) are for the shallow-marine Tarbert formation, the bottom 100 ft (50 layers) are for the fluvial Ness formation. The SPE10 model is structurally simple but highly heterogeneous in porosity and permeability Fig. 3. It poses significant challenges to numerical simulators.

Fig. 4.

Example 3 (SPE10): Numerical pressure for coarse, medium, and fine meshes

The SPE 10 dataset is publicly available at http://www.spe.org/web/csp/. The original data assume the z-axis pointing downwards but use a right-hand coordinate system. A conversion of ordering in blocks is needed for the original data items. We use the code in Matlab Reservoir Simulation Toolbox (MRST) [12] to acquire the needed data.

In this paper, we focus on the Darcy flow part. We use the original permeability data and consider a flow problem. Dirichlet boundary conditions are posed on two boundary faces: $p=1$ for $y=0$; and $p=0$ for $y=1200$. All other four boundary faces have a homogeneous Neumann (no-flow) boundary condition.

We test the WG solver on three meshes (coarse, medium, fine). For better visualization, we tripled the z dimension.

• (i) A coarse mesh with $12\,\times \,44\,\times \,17$ partitions. For the WG$(Q_0,Q_0;RT_{[0]})$ solver, there are 8, 976 pressure degrees of freedom (DOFs) for element interiors; 28, 408 pressure DOFs for all faces, and 37, 384 total DOFs. The local $RT_{[0]}$ spaces are used to compute the discrete weak gradients of the pressure basis functions, but they do not constitute any DOFs.

• (ii) A medium mesh with $30 \times 110 \times 85$ partitions. We use WG$(Q_0,Q_0;RT_{[0]})$ again. There are 280, 500 DOFs for the pressure in element interiors, 856, 700 DOFs for all faces, and totally 1, 137, 200 (about 1M) DOFs.

• (iii) A fine mesh with $60 \times 220 \times 85$ partitions, which is the same as the original gridblock. Again WG$(Q_0,Q_0;RT_{[0]})$ is used. There are 1, 122, 000 interior DOFs; 3, 403, 000 face DOFs; and a total 4, 525, 000 (about 4M) DOFs.

As shown in Fig. 4, the coarse mesh is too coarse to reveal the reservoir geological features. The medium-mesh result is good enough to reflect the channel features of the fluvial Ness formation. The fine-mesh result is smoother and exposes further details about the heterogeneity, but requires expensive computation. Tables 3 and 4 (results on a server with 40 Intel CPUs) together demonstrate that our new WG solver is more efficient than the classical MFEM.

Table 3. Example 3: SPE10_Darcy by WG$(Q_0,Q_0;RT_{[0]})$ on 3 meshes

Table 4. Example 3: SPE10_Darcy by MFEM$(RT_{[0]},Q_0)$ on 3 meshes

6 Concluding Remarks

The novel weak Galerkin finite element methods represent a different type of methodology for solving many real-world problems modeled by partial differential equations. There have been efforts on implementing WG FEMs in Matlab and C++. But the work reported in this paper represents the first ever attempt for implementing WG FEMs in a popular finite element package like deal.II. This shall provide open access to the scientific community for examining usefulness of the WG methodology for large-scale scientific computing tasks.

Listed below are some projects for further research.

• (i) Preconditioning and parallelization of the WG solver for Darcy flow;

• (ii) deal.II implementation for coupled WG Darcy solvers and transport solvers for the full problem of SPE10 and alike;

• (iii) deal.II implementation for both 2-dim and 3-dim for the 2-field poroelasticity solver developed in [10];

• (iv) Implementation of WGFEMs for triangular/tetrahedral meshes on FEniCS;

• (v) Comparison with the hybridizable discontinuous Galerkin (HDG) methods.