A discontinuous Galerkin discretization of elliptic problems with improved convergence properties using summation by parts operators

TL;DR

This paper introduces an efficient discontinuous Galerkin discretization for elliptic problems that leverages summation by parts operators to achieve improved convergence properties, especially for gradients, blending hyperbolic and elliptic perspectives.

Contribution

It demonstrates how to implement a DG method with SBP operators for elliptic problems, enhancing convergence and computational efficiency.

Findings

High order convergence of gradients achieved

Efficient implementation using SBP operators

Improved error estimates for elliptic solutions

Abstract

Nishikawa (2007) proposed to reformulate the classical Poisson equation as a steady state problem for a linear hyperbolic system. This results in optimal error estimates for both the solution of the elliptic equation and its gradient. However, it prevents the application of well-known solvers for elliptic problems. We show connections to a discontinuous Galerkin (DG) method analyzed by Cockburn, Guzm\'an, and Wang (2009) that is very difficult to implement in general. Next, we demonstrate how this method can be implemented efficiently using summation by parts (SBP) operators, in particular in the context of SBP DG methods such as the DG spectral element method (DGSEM). The resulting scheme combines nice properties of both the hyperbolic and the elliptic point of view, in particular a high order of convergence of the gradients, which is one order higher than what one would usually expect from DG methods for elliptic problems.

Tables

Table 1. Table 1 : Summary of the numerical experiment setups.

Setup	Dim.	Domain	Solution	Boundary Condition
Setup 1	1D	$(-1,1)$	$\phi \left(x\right)=\exp \left(-10x^2\right)$	Dirichlet
Setup 2	1D	$(-2,2)$	$\phi \left(x\right)=\exp \left(-10x^2\right)-\sqrt{\pi /10}\mathrm{erf}\left(2\sqrt{10}\right)/4$	periodic
Setup 3	2D	${(-0.5,0.5)}^2$	$\phi (x,y)=\cos \left(\pi x\right)\cos \left(\pi y\right)$	Dirichlet
Setup 4	2D	${(-1,1)}^2$	$\phi (x,y)=2\cos \left(\pi x\right)\sin \left(2\pi y\right)$	periodic

Table 2. Table 2: Numerical results of the 1D convergence experiments summarized in Table  1 .

$N$	Error $\phi$	EOC	Error $q$	EOC
10	2.42e-02		6.88e-02
20	3.16e-03	2.94	8.64e-03	2.99
40	3.97e-04	2.99	1.08e-03	3.01
80	4.96e-05	3.00	1.34e-04	3.00
160	6.19e-06	3.00	1.67e-05	3.00

Table 3. Table 3: Numerical results of the 2D convergence experiments summarized in Table  1 .

$N$	Error $\phi$	EOC	Error $q_1$	EOC	Error $q_2$	EOC
4	7.49e-02		3.34e-01		3.34e-01
8	2.51e-03	4.90	1.15e-02	4.86	1.15e-02	4.86
16	1.70e-04	3.88	8.48e-04	3.77	8.48e-04	3.77
32	1.46e-05	3.54	8.26e-05	3.36	8.26e-05	3.36
64	1.46e-06	3.33	9.23e-06	3.16	9.23e-06	3.16

Equations

$${\partial }_t\phi -\nabla \cdot q$$

$${\partial }_{t}q-\frac{1}{T_r}\nabla \phi$$

$${\int }_{x_i}^{x_{i+1}}q{\partial }_x\chi dx={\int }_{x_i}^{x_{i+1}}f\chi dx+{\left[\hat{q}\chi \right]}_{x_i}^{x_{i+1}},{\int }_{x_i}^{x_{i+1}}q\psi dx=-{\int }_{x_i}^{x_{i+1}}\phi {\partial }_x\psi dx+{\left[\hat{\phi }\psi \right]}_{x_i}^{x_{i+1}},$$

$$\hat{\phi }=\{\{\phi \}\}+\frac{\sqrt{T_r}}{2}[[q]],\hat{q}=\{\{q\}\}+\frac{1}{2\sqrt{T_r}}[[\phi ]],$$

$$\hat{\phi }=\{\{\phi \}\}-C_{12}[[\phi ]]+C_{22}[[q]],\hat{q}=\{\{q\}\}+C_{11}[[\phi ]]+C_{12}[[q]],$$

$$q_e$$

$$-D_e{q}_e$$

$$q_l$$

$$q_r$$

$$t_{r,L}^T{q}_r-t_{l,R}^T{q}_l$$

$$-t_{l,R}^T{D}_l{\phi }_l-t_{l,R}^T{M}_l^{-1}{t}_{l,R}\left(\frac{1}{2}(t_{r,L}^T{\phi }_r-t_{l,R}^T{\phi }_l)+\frac{\sqrt{T_r}}{2}(t_{r,L}^T{q}_r-t_{l,R}^T{q}_l)\right).$$

$$t_{r,L}^T{q}_r-t_{l,R}^T{q}_l=-c_1(t_{r,L}^T{\phi }_r-t_{l,R}^T{\phi }_l)+c_2(t_{r,L}^T{D}_r{\phi }_r-t_{l,R}^T{D}_l{\phi }_l),$$

$$c_1=\frac{1}{2}\frac{t_{r,L}^T{M}_r^{-1}{t}_{r,L}-t_{l,R}^T{M}_l^{-1}{t}_{l,R}}{1+\frac{\sqrt{T_r}}{2}(t_{r,L}^T{M}_r^{-1}{t}_{r,L}+t_{l,R}^T{M}_l^{-1}{t}_{l,R})},c_2=\frac{1}{1+\frac{\sqrt{T_r}}{2}(t_{r,L}^T{M}_r^{-1}{t}_{r,L}+t_{l,R}^T{M}_l^{-1}{t}_{l,R})}.$$

$$q_l$$

$$q_r$$

$$q_e=D_e{\phi }_e$$

$$+M_e^{-1}{t}_{e,L}\left(\frac{1}{2}(1+\sqrt{T_r}{c}_{1,L})[[\phi {]]}_L-\frac{\sqrt{T_r}}{2}{c}_{2,L}[[D\phi {]]}_L\right).$$

$$\hat{\phi }={\phi }^{\mathrm{boundary}},\hat{q}=q^{\mathrm{interior}}+\frac{1}{2\sqrt{T_r}}[[\phi ]].$$

$$q_e=D_e{\phi }_e+M_e^{-1}{t}_{e,R}\left(\frac{1}{2}(1-\sqrt{T_r}{c}_{1,R})[[\phi {]]}_R+\frac{\sqrt{T_r}}{2}{c}_{2,R}[[D\phi {]]}_R\right)+M_e^{-1}{t}_{e,L}(t_{e,L}^T{\phi }_e-{\phi }^{\mathrm{BC}}).$$

$$q_e=D_e{\phi }_e+M_e^{-1}{t}_{e,R}({\phi }^{\mathrm{BC}}-t_{e,R}^T{\phi }_e)+M_e^{-1}{t}_{e,L}\left(\frac{1}{2}(1+\sqrt{T_r}{c}_{1,L})[[\phi {]]}_L-\frac{\sqrt{T_r}}{2}{c}_{2,L}[[D\phi {]]}_L\right).$$

$$L_r=\frac{1}{2\pi }\frac{(x_{\max }-x_{\min })(y_{\max }-y_{\min })}{\sqrt{(x_{\max }-x_{\min }{)}^2+(y_{\max }-y_{\min }{)}^2}}$$

Full text

**AMS subject classification. 65N30, 65N35, 65N06, 65M60, 65M70, 65M06 **

1 Introduction

Solving a Poisson equation $-\Delta\varphi=f$ in a bounded domain $\Omega\subset\mathbb{R}^{d}$ with appropriate boundary conditions (BCs) is a key task in many scientific simulations. Nishikawa [17] proposed to compute numerical solutions as steady state limits of the hyperbolic system

[TABLE]

where $T_{r}>0$ is a relaxation time that can be chosen to accelerate the convergence to the steady state [19]. Some earlier works on this “hyperbolic heat equation” are [6, 16, 13]; some later articles based on the idea are [18, 2, 1, 8].

The “hyperbolic diffusion” approach enables optimal convergence not only of the potential $\varphi$ but also of the gradient $q$ for discontinuous Galerkin (DG) methods. Moreover, it simplifies the coupling to hyperbolic equations in multi-physics problems such as astrophysical fluid flows with self-gravity [24]. However, it would be nice to keep the superconvergence properties of the gradients in a classical elliptic formulation to use state-of-the-art high-performance solvers [12, 9]. Focusing on DG methods, we will explain that the steady-state formulation of (1.1) is equivalent to a scheme analyzed by [7]. However, this method appears to be difficult to implement since a linear system needs to be solved to compute the gradient $q$, i.e., to evaluate the residual of the elliptic discretization. We will explain how this difficulty can be solved for methods using summation by parts (SBP) operators, in particular for discontinuous Galerkin spectral element methods (DGSEM) using Gauss-Lobatto-Legendre nodes [11]. See [14] for more observations how the SBP structure of DGSEM can be used to analyze and improve DG methods for elliptic problems.

2 Main result

We focus on 1D for simplicity. All results extend to the multi-dimensional case using tensor product spaces, e.g., DGSEM. The weak formulation of the steady state of (1.1) with test functions $\psi,\chi$ on an interval $(x_{i},x_{i+1})$ is

[TABLE]

where $\widehat{q},\widehat{\varphi}$ are numerical fluxes. This steady state formulation is obtained by multiplying the second equation by $T_{r}$. Thus, the system depends on the relaxation time $T_{r}$ only via the numerical fluxes $\widehat{q},\widehat{\varphi}$ associated to the time-dependent problem (1.1). The classical upwind numerical fluxes of the hyperbolic system (1.1) are [19]

[TABLE]

where ${\{\mkern-6.0mu\{}\cdot{\}\mkern-6.0mu\}}$ denotes the arithmetic mean and ${[\mkern-3.0mu[}\cdot{]\mkern-3.0mu]}$ the jump at an interface. In the context of the Poisson equation $-\Delta\varphi=f$, these numerical fluxes fit into the classical framework of [3]; they are consistent and conservative (and hence result in adjoint consistency). [5] analyzed similar numerical fluxes of the form

[TABLE]

which match (2.2) for $C_{12}=0$. Using polynomials of degree $p$, they proved that $(\varphi,q)$ converge with orders $(p+1,p+1/2)$ on general meshes if $|C_{12}|$, $|C_{11}|$, and $C_{22}$ are of order unity. [7] extended this analysis and proved that $(\varphi,q)$ converge with optimal order $p+1$ if $C_{11}>0$, $C_{22}>0$, $C_{11}\propto 1/C_{22}$, and $C_{11},|C_{12}|$ are bounded, which is exactly the situation for the numerical fluxes (2.2). However, [7] noted “Of course, the DG methods under consideration are difficult to implement” since the numerical flux $\widehat{\varphi}$ used to compute the gradient $q$ depends on $q$. Thus, a linear system needs to be solved to compute the residual of the elliptic discretization.

In each element $e$, SBP operators [25, 10] are given by i) a discrete derivative operator $D_{e}$ approximating $\partial_{x}$, ii) a diagonal mass/norm matrix $M_{e}$ approximating the $L^{2}$ inner product, and iii) interpolation operators $t_{e,L/R}^{T}$ evaluating a numerical solution at the left/right boundary of the element $e$. Throughout, we use nodal approximations with grid nodes at the boundaries of each element as in DGSEM. In this case, $t_{e,L}^{T}=(1,0,\dots,0)$ and $t_{e,R}^{T}=(0,\dots,0,1)$. Furthermore, we require the SBP condition $M_{e}D_{e}+(M_{e}D_{e})^{T}=t_{e,R}t_{e,R}^{T}-t_{e,L}t_{e,L}^{T}$ mimicking integration by parts [11].

In this framework, the DG discretization (2.1) can be written in the equivalent strong form as

[TABLE]

Here, $q_{e}$, $\varphi_{e}$, and $f$ are the vectors of coefficients representing the respective polynomials in element $e$ in the chosen nodal basis. Now, we are prepared to formulate the main result of this short note.

Theorem 2.1**.**

The discrete gradient $q$ of (2.1) with numerical fluxes (2.2) can be evaluated locally using only surface values of $\varphi$ and $D\varphi$ from neighboring elements if diagonal-norm SBP operators including the boundaries are used. In particular, DGSEM is included in this class of SBP operators.

Proof.

It suffices to consider two elements indicated by subscripts $l,r$. Abbreviating surface terms not belonging to their common interface as $ST_{l,L}$ (left element, left surface) and $ST_{r,R}$ (right element, right surface), the corresponding discretizations are

[TABLE]

Again, $q_{l/r}$ and $\varphi_{l/r}$ are the vectors of coefficients representing the respective polynomials in the left/right element in the chosen nodal basis. Since boundary nodes are included and the mass matrix is diagonal, the surface terms vanish everywhere except at their corresponding interface nodes. In particular, the restriction of $q_{l}$ to the right surface of element $l$ is not influenced by the left surface term $ST_{l,L}$. Hence, the jump of boundary values of $q$ is

[TABLE]

This equation can be solved for the jump of interface values of $q$,

[TABLE]

where

[TABLE]

Note that $c_{1}$ vanishes for uniform grids with symmetric quadrature rules. In general, $c_{1}$ and $c_{2}$ depend on the grid spacing. Inserting this expression of the jump of $q$ at the interface into (2.5) yields

[TABLE]

To sum up, the gradient in an interior element $e$ can be computed explicitly as

[TABLE]

At a boundary point where a Dirichlet condition is imposed weakly for the potential $\varphi$, the numerical fluxes based on an energy analysis for the hyperbolic diffusion system are the ones used by [7], i.e.,

[TABLE]

Hence, the numerical flux $\widehat{\varphi}$ does not depend on $q$ at a Dirichlet boundary and no special care is needed. Thus, the discretization in the element at the left boundary is

[TABLE]

Similarly, the gradient in the element at the right boundary is

[TABLE]

Thus, the gradient can be computed locally using surface values of $\varphi$ and $D\varphi$ from neighboring elements. ∎

3 Numerical experiments

We demonstrate the convergence properties of the method for several Poisson problems summarized in Table 1. The right-hand side $f$ is chosen based on the solution $\varphi$. We use Dirichlet BCs for non-periodic setups and vanishing mean values of $\varphi$ for periodic BCs. The non-periodic 2D setup is taken from [7]. We choose the relaxation time $T_{r}$ as recommended in [19], i.e., $T_{r}=L_{r}^{2}$ where the reference length scale $L_{r}$ is set to $L_{r}=(x_{\mathrm{max}}-x_{\mathrm{min}})/(2\pi)$ for a 1D interval $(x_{\mathrm{min}},x_{\mathrm{max}})$ and

[TABLE]

for a 2D rectangle $(x_{\mathrm{min}},x_{\mathrm{max}})\times(y_{\mathrm{min}},y_{\mathrm{max}})$.

All methods are implemented in Julia [4]. We use Trixi.jl [23, 24] to compute steady state solutions of the hyperbolic system (1.1) and SummationByPartsOperators.jl [22] to implement the corresponding elliptic approach. All source code required to reproduce the numerical experiments is available online [21].

We compute the $L^{2}$ errors of the numerical solutions on meshes with $N$ elements per coordinate direction using the Gauss-Lobatto-Legendre quadrature rule associated with the DGSEM operators. The 1D results including the experimental order of convergence (EOC) are shown in Table 2. These results are computed using direct sparse solvers distributed with Julia [4]. Clearly, both the potential $\varphi$ and the gradient $q$ converge with optimal order $p+1$ for polynomials of degree $p$. The results obtained by evolving the hyperbolic system (1.1) match the results obtained from the elliptic implementation and are thus not shown.

The 2D results are shown in Table 3. These results are computed using the conjugate gradients (CG) implementation of Krylov.jl [15] using matrix-free operators based on the interface of LinearOperators.jl [20]. Again, both the potential $\varphi$ and the gradient $q$ converge with optimal order.

Acknowledgments

Special thanks to Jesse Chan for discussions related to this manuscript and comments on an early draft.