A discrete maximum principle for the weak Galerkin finite element method on nonuniform rectangular partitions

In this paper, we are concerned with the development of a discrete maximum principle for the weak Galerkin finite element approximations of convection-diffusion-reaction problems. For simplicity, we consider the problem of seeking an unknown function $u=u(x)$ satisfying

where $\Omega$ is a bounded polytopal domain in $\mathbb{R}^{d}\;(d\geq 2)$ with boundary $\partial\Omega$, $\alpha=\alpha(x)$ is the diffusion coefficient, ${\bm{\beta}}={\bm{\beta}}(x)$ is the convection, and $c=c(x)$ is the reaction coefficient in relevant applications. We assume that $\alpha$ is sufficiently smooth, ${\bm{\beta}}\in[W^{1,\infty}(\Omega)]^{d}$, and $c$ is piecewise constant with respect to a partition of the domain. For well-posedness of the problem (1)-(2), we assume $f=f(x)\in L^{2}(\Omega)$, $g=g(x)\in H^{\frac{1}{2}}(\partial\Omega)$, and

for a constant $\alpha_{0}>0$.

The weak Galerkin (WG) finite element method, recently introduced in [20, 21, 15], refers to a natural extension of the standard finite element methods [3, 8] where the differential operators are approximated as discrete distributions or discrete weak derivatives and the element continuity circumvented by properly selected stabilizers. The method has good flexibility in making use of discontinuous elements while sharing the simple formulation of the classical continuous or conforming finite element methods. Weak Galerkin can be easily implemented on finite element partitions consisting of general polygonal or polyhedral elements through the usual assembling strategy of element stiffness matrices. The method has gained a lot attention and popularity recently and has been successfully applied to a variety of partial differential equations, see [20, 11, 14, 19] and the references therein for more details.

A simplified formulation of the weak Galerkin finite element method has been developed in [10, 13] for second order elliptic equations in conjunction with the study of superconvergence and error estimates. This simplification idea was further applied to the Stokes equation in [12] for the development of a superconvergence theory on nonuniform rectangular partitions for both the velocity and the pressure approximations. In the simplified weak Galerkin (SWG), the degrees of freedom associated with the unknowns in the interior of each element are eliminated from the usual weak Galerkin method, yielding a numerical scheme with significantly reduced computational complexity. For pure diffusion equations (i.e., ${\bm{\beta}}=0$ and $c=0$ in (1)), the simplified weak Galerkin is equivalent to the usual weak Galerkin in the sense that the numerical approximations on the element boundary are the same. But for the full convection-diffusion-reaction equation, these two methods give different numerical solutions on the element boundary. Like WG, the simplified weak Galerkin preserves the important mass conservation property locally on each element and allows the use of general polygonal partitions.

The convection-diffusion-reaction equation (1) is known to satisfy the following maximum principle: If $u\in C^{2}(\Omega)\cap C^{1}(\bar{\Omega})$ is the solution of (1)-(2) with $\alpha,{\bm{\beta}},c\in C^{1}(\Omega)$ and $f\leq 0$ in $\Omega$, then $u$ attains its maximum value on $\partial\Omega$ if $c=0$ and $u$ attains its non-negative maximum value on $\partial\Omega$ if $c\geq 0$ [7]. A discrete maximum principle (DMP) refers to a similar statement for numerical solutions of (1)-(2) on specific grid points. The discrete maximum principle is of great importance from physical point of views in scientific computing (e.g. helping avoid non-physical numerical solutions). In the last two decades, a great deal of effort has been devoted to the search of DMP-preserving numerical methods, see [1, 2, 5, 6, 16, 18]. For isotropic diffusion problems, it was shown by Ciarlet and Raviart [4] that the $P_{1}$-conforming finite element solutions satisfy a DMP if all the triangular elements have non-obtuse dihedral angles. This nonobtuse-angle condition was improved in [17] by a weaker condition (namely, the Delaunay condition) that requires the sum of any pair of angles facing a common interior edge be less than or equal to $\pi$. A more recent work on DMP was developed in [23] for $P_{1}$-conforming finite element approximations of quasi-linear second order elliptic equations by using the de Giorgi approach developed in PDE analysis. For discontinuous Galerkin finite element approximations, some DMP-preserving schemes were developed in [24] for convection-diffusion equations on triangular meshes. In the weak Galerkin context and for anisotropic diffusion problems, a DMP result was established for the lowest order WG solutions without stabilization in [9]. In [22], the authors developed a DMP theory for the WG numerical solutions on the element boundary under a weak acute angle condition for triangular partitions. For weak Galerkin on rectangular partitions, a DMP result was reported numerically in [22], but no theory was developed over there.

The goal of this paper is to establish a discrete maximum principle for the simplified weak Galerkin finite element approximations of (1)-(2) on nonuniform rectangular partitions. We show that a DMP is satisfied by the numerical solutions arising from SWG with certain values of the stabilization parameter on rectangular partitions for elements with aspect ratio in $[0.5,2]$. For a better understanding of the SWG, we shall compute its global stiffness matrix on uniform rectangular partitions, and develop a 5- and 7-point finite difference scheme for the model diffusion equation. As a variant of the weak Galerkin scheme, these finite difference methods can be shown to preserve the important mass conservation property locally on each cell. We note that a similar finite difference method has been developed in [12]) for the Stokes equation.

The paper is organized as follows: In Section 2, we state the simplified weak Galerkin finite element method for the model problem (1)-(2). In Section 3, we establish a comprehensive DMP theory for the SWG approximations. Section 4 is devoted to the derivation of a technical inequality useful to the DMP development. In Section 5, we devise a finite difference scheme for the diffusion equation on uniform Cartesian grids based on the SWG formulation. Finally, in Section 6, we report some numerical results for a verification of the discrete maximum principle.

In the rest of the paper, we assume $d=2$ and shall use the standard notations for Sobolev spaces and norms [3, 8]. For any open set $D\subset\mathbb{R}^{2}$, $\|\cdot\|_{s,D}$ and $(\cdot,\cdot)_{s,D}$ denote the norm and inner-product in the Sobolev space $H^{s}(D)$ consisting of square integrable partial derivatives up to order $s$. When $s=0$ or $D=\Omega$, we shall drop the corresponding subscripts in the norm and inner-product notation.

Let ${\mathcal{T}}_{h}=\{T\}$ be a shape-regular polygonal partition of the domain $\Omega$. For $T\in{\mathcal{T}}_{h}$, denote by $h_{T}$ its diameter and by $N$ the number of edges. The meshsize of ${\mathcal{T}}_{h}$ is defined as $h=\max_{T\in{\mathcal{T}}_{h}}h_{T}$. For each edge $e_{i},\ i=1,\ldots,N$, let $M_{i}$ be the midpoints and ${\bf n}_{i}$ be the outward normal direction of $e_{i}$; see Fig. 1.

Let $v_{b}$ be a piecewise constant function defined on the boundary of $T$. The weak gradient of $v_{b}$ [21, 15] is given by

where $v_{b,i}=v_{b}|_{e_{i}}$, $|e_{i}|$ is the length of the edge $e_{i}$, and $|T|$ is the area of $T$. The weak gradient $\nabla_{w}v_{b}$ can be seen to satisfy the following equation:

for all constant vector $\bm{\phi}$, where and in what follows of the paper, $\langle\cdot,\cdot\rangle_{\partial T}$ is the notation for the usual inner product in $L^{2}({\partial T})$.

Denote by $W(T)$ the local finite element space consisting of piecewise constant functions on ${\partial T}$. The global finite element space, denoted by $W({\mathcal{T}}_{h})$, is given by patching all the local elements $W(T)$ through common values on interior edges. Denote by $W_{h}^{0}({\mathcal{T}}_{h})$ the closed subspace of $W_{h}({\mathcal{T}}_{h})$ consisting of functions with vanishing boundary values.

Let ${P}_{k}(T)$ be the space of polynomials of degree $k\geq 0$ on $T$. Each $v_{b}\in W(T)$ can be extended to $T$ as a linear function ${\mathfrak{s}}(v_{b})\in{P}_{1}(T)$ by the following equations:

It is not hard to show that ${\mathfrak{s}}(u_{b})$ is well defined by (6).

On each element $T\in{\mathcal{T}}_{h}$, we introduce three bilinear forms:

Furthermore, let

for $u_{b},v_{b}\in W(T)$. To enforce a weak continuity, we introduce the following stabilizer:

where $Q_{b}$ is the $L^{2}$ projection operator onto $W(T)$. It is clear that $Q_{b}u$ is the average of $u$ on each edge. With an abuse of notation, but without confusion, we use $Q_{b}(g)$ to denote the average of the Dirichlet value $g$ on each boundary edge.

The following result on the solution existence and uniqueness has been established by the authors in [13].

The model problem (1)-(2) is known to satisfy the following maximum principle: If $u\in C^{2}(\Omega)\cap C^{1}(\bar{\Omega})$ is the solution of (1)-(2) with $\alpha,{\bm{\beta}},c\in C^{1}(\Omega)$ and non-positive $f\in C(\Omega)$, then $u$ attains its maximum value on $\partial\Omega$ when $c=0$ and attains its non-negative maximum value on $\partial\Omega$ when $c\geq 0$ [7]. In this section, we show that the maximum principle also holds true for the numerical solutions arising from the SWG scheme (12) on nonuniform rectangular partitions. From now on, the finite element partitions ${\mathcal{T}}_{h}$ is assumed to contain only rectangular elements.

In practical computation, the load linear form $\sum_{T}(f,{\mathfrak{s}}(v_{b}))_{T}$ in the SWG scheme (12) must be approximated by using numerical integrations on each element $T$. Let us first derive a discrete version for $\sum_{T}(f,{\mathfrak{s}}(v_{b}))_{T}$ on which the maximum principles will be established.

Let $T\in{\mathcal{T}}_{h}$ be a rectangular element depicted in Fig. 2 with center $M_{T}=(x_{T},y_{T})$. Denote by $h_{x}:=|e_{3}|=|e_{4}|$ and $h_{y}:=|e_{1}|=|e_{2}|$ the local meshsize in $x$ and $y$ directions. From (6), the linear extension of $u_{b}\in W(T)$ can be represented as (see Lemma 6.1 in [10] for details)

where

It follows that

Let $\phi_{i}$ be the local basis function associated with the edge $e_{i}$ (i.e., $\phi_{i}=1$ on $e_{i}$ and $\phi_{0}=0$ on other edges). For any $v_{b}=\sum_{i=1}^{4}v_{b,i}\phi_{i}\in W(T)$, we have from (15)

where

It follows that

Note that ${\mathfrak{s}}(\phi_{1})$ is linear in $x$- direction, and constant in $y$-. Thus, the integral $(f,{\mathfrak{s}}(\phi_{1}))_{T}$ can be approximated by using the mid-point rule in $y$-direction and the Simpson’s rule in $x$-direction:

From the Taylor expansion, we have

Substituting the above into (21) yields

where $\sigma=h_{x}/h_{y}$. Analogously, we have

Using (20)-(25), we may approximate $(f,{\mathfrak{s}}(v_{b}))_{T}$ by $(f,{\mathfrak{s}}(v_{b}))_{h,T}$ given as follows:

In the case that $f=f(x,y)$ is non-smooth, the values $f(M_{i})$ and $f(x_{T},y_{T})$ can be substituted by a local average of $f$ around $M_{i}$ and $(x_{T},y_{T})$, respectively. Our computational version of the SWG finite element method for the elliptic equation (1)-(2) can be stated as follows:

We are now in a position to state the main result of this section on discrete maximum principles.

The goal of this section is to verify the validity of the assumption (32). From $\theta=(u_{b}-k)^{+}$ and $\psi=(u_{b}-k)^{-}$, we may rewrite the assumption as follows:

Let $T\in{\mathcal{T}}_{h}$ be a rectangular element depicted in Fig. 2. For any $v_{b}\in W(T)$, define $v_{b}^{+}$ and $v_{b}^{-}$ as follows

It is easy to see that

The following lemma provides a version of (34) on the element $T$.