Manycore parallel computing for a hybridizable discontinuous Galerkin nested multigrid method

TL;DR

This paper develops a parallel hybridizable discontinuous Galerkin nested multigrid solver optimized for many-core processors, achieving high efficiency and convergence rates for high polynomial orders in computational PDEs.

Contribution

It introduces a fine-grain parallelization strategy for HDG nested multigrid methods on many-core architectures, utilizing matrix-free techniques and data locality for high performance.

Findings

Achieves ideal convergence rates of 0.2 or less for high polynomial orders.

Attains 80% of peak bandwidth performance with high-order polynomials.

Demonstrates high efficiency and scalability on Intel Xeon Phi (Knights Landing) processors.

Abstract

We present a parallel computing strategy for a hybridizable discontinuous Galerkin (HDG) nested geometric multigrid (GMG) solver. Parallel GMG solvers require a combination of coarse-grain and fine-grain parallelism to improve time to solution performance. In this work we focus on fine-grain parallelism. We use Intel's second generation Xeon Phi (Knights Landing) many-core processor. The GMG method achieves ideal convergence rates of $0.2$ or less, for high polynomial orders. A matrix free (assembly free) technique is exploited to save considerable memory usage and increase arithmetic intensity. HDG enables static condensation, and due to the discontinuous nature of the discretization, we developed a matrix vector multiply routine that does not require any costly synchronizations or barriers. Our algorithm is able to attain 80\% of peak bandwidth performance for higher order polynomials. This is possible due to the data locality inherent in the HDG method. Very high performance is realized for high order schemes, due to good arithmetic intensity, which declines as the order is reduced.

Tables

Table 1. Table 1: Errors and convergence rates of the HDG scheme, on a Cartesian mesh of N × N 𝑁 𝑁 N\times N elements.

		${\Vert u_h-u\Vert }_{L^2\left(\Omega \right)}$		${\Vert u_h^{\ast }-u\Vert }_{L^2\left(\Omega \right)}$		${\Vert q_h-q\Vert }_{L^2\left(\Omega \right)}$		${\Vert q_h^{\ast }-q\Vert }_{L^2\left(\Omega \right)}$
$p$	$N$	Error	Rate	Error	Rate	Error	Rate	Error	Rate
6	2	1.05e-07	-	3.02e-09	-	2.39e-07	-	1.12e-07	-
	4	8.51e-10	6.95	1.13e-11	8.06	1.95e-09	6.94	8.41e-10	7.06
	8	6.97e-12	6.93	4.45e-14	7.99	1.61e-11	6.92	6.65e-12	6.98
	16	5.60e-14	6.96	7.12e-16	5.97	1.41e-13	6.84	6.36e-14	6.71
	32	8.79e-16	5.99	5.35e-16	4.12	2.11e-13	-5.82	1.47e-13	-1.21
5	2	7.53e-07	-	3.50e-08	-	1.77e-06	-	9.80e-07	-
	4	1.69e-08	5.47	3.03e-10	6.85	3.91e-08	5.50	1.81e-08	5.76
	8	2.90e-10	5.87	2.43e-12	6.96	6.72e-10	5.86	2.95e-10	5.94
	16	4.72e-12	5.94	1.92e-14	6.98	1.10e-11	5.94	4.69e-12	5.97
	32	7.52e-14	5.97	3.68e-16	5.70	2.05e-13	5.74	1.01e-13	5.53
4	2	1.46e-05	-	7.04e-07	-	3.32e-05	-	1.79e-05	-
	4	4.92e-07	4.89	1.09e-08	6.02	1.13e-06	4.88	5.49e-07	5.03
	8	1.65e-08	4.90	1.73e-10	5.97	3.80e-08	4.89	1.76e-08	4.96
	16	5.37e-10	4.94	2.75e-12	5.98	1.24e-09	4.94	5.62e-10	4.97
	32	1.71e-11	4.97	4.33e-14	5.99	3.98e-11	4.97	1.78e-11	4.98
3	2	8.61e-05	-	7.58e-06	-	2.05e-04	-	1.35e-04	-
	4	7.78e-06	3.47	2.63e-07	4.85	1.80e-05	3.52	9.86e-06	3.78
	8	5.49e-07	3.82	8.63e-09	4.93	1.27e-06	3.82	6.59e-07	3.90
	16	3.63e-08	3.92	2.77e-10	4.96	8.43e-08	3.91	4.28e-08	3.95
	32	2.33e-09	3.96	8.80e-12	4.98	5.43e-09	3.96	2.73e-09	3.97

Table 2. Table 2: KNL testbed specifications.

Processor number	7210
# of cores	64
Processor base frequency	1.30 GHz
Cache	32 MB L2
RAM	384 GB DDR4
MCDRAM	16 GB
Instruction set	64-bit
Instruction set extension	Intel® AVX-512
Operating system	CentOS 7.2
Compiler	Intel® Parallel Studio XE 16

Table 3. Table 3: Arithmetic intensity and GFLOP/s for the local solvers.

Polynomial Order	AI	GFLOP/s
0	0.84	0.0032
1	1.41	1.15
2	2.60	5.129
3	4.33	21.18
4	6.56	53.66

Equations

$$-\nabla \cdot (\mathbf{K}\nabla u)$$

$$u$$

$$\mathbf{q}$$

$${\mathbf{K}}^{-1}\mathbf{q}+\nabla u$$

$$\nabla \cdot \mathbf{q}$$

$$u$$

$$(u,v{)}_{{\mathcal{E}}_h}$$

$$⟨\mathbf{w}\cdot \mathbf{n},v{⟩}_{{\Gamma }_h}$$

$$W_h$$

$$M_h$$

$$({\mathbf{q}}_h,\mathbf{v}{)}_{{\mathcal{E}}_h}-(u_h,\nabla \cdot \mathbf{v}{)}_{{\mathcal{E}}_h}+⟨{\lambda }_h,\mathbf{v}\cdot \mathbf{n}{⟩}_{{\Gamma }_h^o}$$

$$-({\mathbf{q}}_h,\nabla w{)}_{{\mathcal{E}}_h}+⟨{\mathbf{q}}_h\cdot \mathbf{n},w{⟩}_{{\Gamma }_h}+⟨\tau (u_h-{\lambda }_h),w{⟩}_{{\Gamma }_h^o}+⟨\tau u_h,w{⟩}_{\partial \Omega }$$

$$⟨{\mathbf{q}}_h\cdot \mathbf{n},\mu {⟩}_{{\Gamma }_h}+⟨\tau (u_h-{\lambda }_h),\mu {⟩}_{{\Gamma }_h^o}+⟨\tau u_h,\mu {⟩}_{\partial \Omega }$$

$${\int }_e{P}_h{g}_D\mu ={\int }_e{g}_D\mu ,\forall \mu \in {\mathbb{Q}}_p(e),\forall e\in \partial \Omega .$$

$${\mathbf{A}}_K{\mathbf{Q}}_K-{\mathbf{B}}_K^T{\mathbf{U}}_K+{\mathbf{C}}_K{\mathbf{\Lambda }}_K$$

$${\mathbf{B}}_K{\mathbf{Q}}_K+{\mathbf{D}}_K{\mathbf{U}}_K+{\mathbf{E}}_K{\mathbf{\Lambda }}_K$$

$$\begin{bmatrix}\bm{Q}_{K}\\ \bm{U}_{K}\end{bmatrix}=\begin{bmatrix}{\bm{A}}_{K}&-{\bm{B}}_{K}^{T}\\ {\bm{B}}_{K}&{\bm{D}}_{K}\end{bmatrix}^{-1}\Bigg{(}-\begin{bmatrix}{\bm{C}}_{K}\\ {\bm{E}}_{K}\end{bmatrix}\bm{\Lambda}_{K}+\begin{bmatrix}\bm{R}_{K}\\ \bm{F}_{K}\end{bmatrix}\Bigg{)}.$$

$${\mathbf{G}}_{K_{1}e}{\mathbf{Q}}_{K_1}+{\mathbf{G}}_{K_{2}e}{\mathbf{Q}}_{K_2}+\tau {\mathbf{H}}_{K_{1}e}{\mathbf{U}}_{K_1}+\tau {\mathbf{H}}_{K_{2}e}{\mathbf{U}}_{K_2}-2\tau {\mathbf{M}}_e{\mathbf{\Lambda }}_{K_1{K}_2}=0.$$

$${\mathbf{G}}_{Ke}{\mathbf{Q}}_K+\tau {\mathbf{H}}_{Ke}{\mathbf{U}}_K={\mathbf{S}}_e.$$

$$v(x)=\frac{\sum _{j=0}^n\frac{W_j}{x-x_j}{f}_j}{\sum _{j=0}^n\frac{W_j}{x-x_j}},$$

$$W_j^{\text{GLL}}=(-1{)}^j\sqrt{w_j^{\text{GL}}}\text{leavevmode}\text{leavevmode}\text{leavevmode}{W}_j^{\text{GL}}=(-1{)}^j\sqrt{x_j^{\text{GL}}{w}_j^{\text{GL}}},$$

$$({\xi }_i,{\eta }_j)\stackrel{def}{=}({\xi }_i^{(p)},{\eta }_j^{(p)}),\text{leavevmode}\text{leavevmode}\text{leavevmode}i,j=0,1,\dots ,p,$$

$${\varphi }_I^e(\xi ,\eta )={\ell }_i(\xi ){\ell }_j(\eta ),$$

$$\displaystyle u(x,y)\bigg{|}_{\Omega^{e}}$$

$${\varphi }_I^e({\xi }_m,{\eta }_n)$$

$$(\xi ,\eta )$$

$${\iint }_{\square }u(\xi ,\eta )d\xi d\eta$$

$$\approx \sum _{i=0}^P\sum _{j=0}^P{w}_i{w}_{j}u({\xi }_i,{\eta }_j),$$

$$\mathcal{J}^{e}(\xi,\eta)=\begin{vmatrix}\frac{\partial x}{\partial\xi}&\frac{\partial x}{\partial\eta}\\ \frac{\partial y}{\partial\xi}&\frac{\partial y}{\partial\eta}\end{vmatrix}=\bigg{|}\frac{\partial x}{\partial\xi}\frac{\partial y}{\partial\eta}-\frac{\partial y}{\partial\xi}\frac{\partial x}{\partial\eta}\bigg{|}.$$

$${\iint }_{{\Omega }^e}u(x,y)dxdy$$

$$\approx \sum _{i=0}^p\sum _{j=0}^p{w}_i{w}_j{u}^e({\xi }_i,{\eta }_j){\mathcal{J}}_{ij}^e,$$

$${\ell }_j(x)=\frac{\frac{W_j^{\text{GL}}}{x-x_j^{\text{GLL}}}}{\sum _{k=0}^p\frac{W_k^{\text{GL}}}{x-x_k^{\text{GLL}}}}.$$

$$P(u_H)=u_H\forall u_H\in {\Omega }^H.$$

$$(R{v}_h,u_H)=(v_h,u_H),\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall v_h\in ({\Omega }^h{)}^{\prime },\text{leavevmode}\text{leavevmode}\text{leavevmode}\forall u_H\in {\Omega }^H.$$

$$M_k=\{v:\Omega \to \mathbb{R}:v\text{ is continuous }v{|}_{\partial \Omega }=0,\text{leavevmode}v{|}_K\in Q^1(K),\text{leavevmode}\forall K\in {\mathcal{E}}_{k+1}\},$$

$$P_{k}v=\{\begin{array}{ll}v,& \text{mbox}ifk<J,\\ {\Pi }_{M_J}(v{|}_{{\Gamma }_h}),& \text{mbox}ifk=J,\end{array}$$

Full text

{AMS}

65M55,65N30,35J15

1 Introduction

Multigrid methods are among the most efficient solvers for linear systems that arise from the discretization of partial differential equations. The effectiveness of this multilevel technique was first observed by Fedorenko ([21]) in 1964, and popularized by Brandt ([11]) in 1977. Traditionally, multigrid methods have been applied to low order finite difference, finite volume, and continuous finite element approximations discretizations ([53], [56], [25], [12]). The discontinuous Galerkin (DG) discretization brings with it many advantages: can handle complex geometries, has access to $hp$-adaptivity, is capable of satisfying local mass balance (ideal for flow problems and hyperbolic PDEs), and highly parallelizable due to the lack of continuity constraints between elements. However, DG methods often have more degrees of freedom than their continuous counterpart. In addition, high order discretizations rapidly increase the condition number of the discretization operator, which poses a challenge for any linear solver ([7]).

Originally multigrid methods were developed at a time when access to parallelism was limited, and iterative methods that had less concurrency but better convergence rates were favored. The multilevel nature of multigrid can cause challenges in its parallelization; load balancing problems occur due to fine grids have ample data to work with, but coarser levels do not. Moreover, traditional multigrid is a multiplicative method, that is, each level must be completely processed before moving to the next. Additive multigrid methods allow for the simultaneous processing of all levels, but the trade off between concurrency and robustness is often not ideal [8].

Multiplicative multigrid is well known to have sequential complexity, for $N$ data points a single cycle costs $\mathcal{O}(N)$ floating point operations. However, the parallel complexities of V- and F-cycles are polylogarithmic. For multiplicative multigrid, a natural heterogeneous computing strategy is to process fine levels up to a threshold on an coprocessor (or accelerator), and have the remaining coarse levels be processed by coarse-grain parallelism. In this paper all computations are done on a single Xeon Phi coprocessor (Knights Landing), as we only focus on the fine-grain parallelism. However, an offloading model is a natural extension of this work.

In this paper we consider a nested multigrid technique for high order hybridizable discontinuous Galerkin methods; which combines $p$-multigrid and $h$-multigrid. The HDG method is capable of static condensation, which significantly reduces the total number of degrees of freedom compared to classical DG methods. This has important consequences for linear solvers like multigrid, since the cost of multigrid grows with the number of degrees of freedom. Moreover, due to the discontinuous nature of HDG methods, we show that we can efficiently leverage the massive parallelism that many-core processors offer. Through roofline performance modeling and tuning, we show that our implementation results in efficient device utilization as well as significant speedups over a serial implementation. An interesting benefit of the Xeon Phi coprocessor is that it allows us to use traditional parallel programming paradigms like OpenMP, MPI, and pthreads. Consequently, we can take advantage of the massive fine-grain parallelism that the Xeon Phi offers by utilizing these traditional parallel paradigms with significantly limited software intrusion.

The remainder of the paper is organized as follows. Section 2 introduces the model PDE. The relevant finite element notations and HDG discretization is described in section 3. In section 4, we explain how we implement the nodal tensor product basis, quadrature, and barycentric interpolation. Numerical experiments verifying the correct convergences rates for the HDG discretization is conducted in section 5. The core components of the multigrid solver we use is discussed in section 6. This includes the description of intergrid transfer operators, relaxation, coarse grid operators, and the standard multigrid cycle. In section 7 we show through computational experiments that our multigrid method obtains very good convergence rates and error reduction properties. The performance of our algorithm in a parallel setting is evaluated in section 8.

2 Model problem

Consider the elliptic problem

[TABLE]

where $\Omega$ is an open domain in $\mathbb{R}^{2}$ and $\partial\Omega$ denotes the boundary of the domain. Dirichlet datum $g_{D}$ is imposed on the boundary. The function $f$ is the prescribed source function and the matrix $\bf K$ is symmetric positive definite with piecewise constant entries. The div-grad operator in (1) appears in several problems in engineering such as multiphase flow in porous media. It will provide insight into how a DG multigrid solver performs on modern architectures. Moreover, it acts as a gateway to construct very efficient numerical methods for time dependent PDEs that require pressure solves or implicit time stepping.

3 Discretization

HDG methods were designed to address the concern that DG schemes generate more degrees of freedom (DOFs) when compared to continuous Galerkin techniques. For standard DG methods, each degree of freedom is coupled with all other degrees of freedom on a neighboring element. By introducing additional unknowns along element interfaces, the HDG method is able to eliminate all degrees of freedom that do not reside on the interfaces. As such, a significantly smaller linear system is generated, and HDG gains much of its efficiency at higher orders ([37]). It turns out that the HDG method also has a number of attractive properties, namely, the capability of efficient implementations, optimal convergence rates in the potential and flux variables, as well as the availability of a postprocessing technique that results in the superconvergence of the potential variable. HDG methods are a subset of DG methods, so they still retain favorable properties, e.g., local mass balance, ease of $hp$-adaptivity, and the discontinuous nature of the solution variables. A thorough analysis of HDG methods can be found in [18], [20], and [17].

A number of works are available on multigrid for DG methods. Interior penalty methods are the most commonly analyzed, for instance, see [14], [13], [23], [4], and [3]. Most of these works are theoretical, and while they are able to prove convergence, the numerical experiments show rates below what is typically expected from GMG in this model setting. For the interior penalty class of DG methods, it was found that specialized smoothers and careful tuning of the stability parameter was required for better convergence results (see [34], [33]). Local Fourier analysis (local mode analysis, [12]) was applied to interior penalty DG methods in [27], [29], and [28]. However, convergence rates were in the range of $0.4$ to $0.6$ for low order discretiztions $(p\leq 2)$. In addition, local Fourier analysis is applied to a two level multigrid scheme, and is used as a heuristic to estimate GMG performance. Further, it is well known that two grid optimality does not always imply V-cycle optimality (see [44]).

It should be noted that the HDG class of discretizations is quite large, similar to that of standard DG methods. In [20], a unified framework is developed (similar to the work of Arnold et al. in [6]) to create a taxonomy of HDG methods. For instance, one can obtain HDG methods by utilizing one of: Raviart–Thomas DG, Brezzi–Douglas–Marini DG, local DG, or interior penalty DG. In this work we have no need to distinguish between the various HDG methods, because we employ the LDG family of hybridizable methods ([17]). As such, since the LDG method is a mixed technique, one needs to reformulate the underlying equation (1) as a first order system by introducing an auxiliary variable ${\bf q}$:

[TABLE]

We now describe the HDG method. Let $\mathcal{E}_{h}$ be a subdivision of $\Omega$, made of quadrilaterals, $K$, of maximum diameter $h$. The unit normal vector outward of $K$ is denoted by $\bm{n}_{K}$. The mesh skeleton is denoted by $\Gamma_{h}$, that is the union of all the edges. We further decompose the mesh skeleton as $\Gamma_{h}=\Gamma_{h}^{o}\cup\Gamma_{h}^{\partial},$ where $\Gamma_{h}^{\partial}$ denotes the set of all edges on the boundary of the domain, and $\Gamma_{h}^{o}$ the set of all interior edges. The broken Sobolev space is represented by $H^{1}(\mathcal{E}_{h})$; it consists of piecewise $H^{1}$ functions on each mesh element. We use the following short-hand notation for $L^{2}$ inner-product on mesh elements and edges:

[TABLE]

[TABLE]

The underlying approximation spaces for the HDG method are as follows:

[TABLE]

where $\mathbb{Q}_{p}(K)$ is the standard finite element space for quadrilaterals. That is, $\mathbb{Q}_{p}(K)$ is the tensor product of polynomials of degree $p$ on each variable. The same definition (6) applies to $\langle{\bm{w}}\cdot{\bm{n}},\mu\rangle_{\Gamma_{h}}$ for functions ${\bm{w}}\in H^{1}(\mathcal{E}_{h})^{2}$ and $\mu\in M_{h}$.

The HDG method seeks an approximation $({\bm{q}}_{h},u_{h},\lambda_{h})\in{\bm{V}}_{h}\times{W}_{h}\times{M}_{h}$ of the exact solution $({\bm{q}}|_{\Omega},u|_{\Omega},u|_{\Gamma_{h}\setminus\partial\Omega_{D}})$ such that

[TABLE]

for all $({\bm{v}},w,\mu)\in{\bm{V}}_{h}\times{W}_{h}\times{M}_{h}$. The factor $\tau$ is a local stabilization term that is piecewise constant on $\mathcal{E}_{h}$. In the above, $P_{h}g_{D}$ is the $L^{2}$–projection of $g_{D}$, defined by:

[TABLE]

Introducing additional unknowns at first glance does not appear to add much benefit. However, the HDG method allows us to eliminate the unknowns ${\bm{q}}_{h}$ and $u_{h}$ using equations (7) and (9) in an element by element manner to arrive at a weak formulation in terms of $\lambda_{h}$ only. Let $\bm{Q}$, $U$ and $\bm{\Lambda}$ denote the vectors of degrees of freedom for the numerical solutions ${\bm{q}}_{h}$, $u_{h}$ and $\lambda_{h}$ respectively. Since the solutions ${\bm{q}}_{h}$ and $u_{h}$ are discontinuous, we can denote by $\bm{Q}_{K}$ and $\bm{U}_{K}$ the part of the vectors $\bm{Q}$, $\bm{U}$ corresponding to the degrees of freedom located on $K$. We also denote by $\bm{\Lambda}_{K}$ the part of the vector $\bm{\Lambda}$ that corresponds to the degrees of freedom located on the boundary of $K$. We can express equations (7), (8) in matrix form:

[TABLE]

Further, one can obtain the local degrees of freedom:

[TABLE]

The above inverse is well defined and elimination of the local degrees of freedom for $\bm{q}$ and $u$ can be done in parallel, independently of one another [20]. Equation (9) can be written in matrix form. Fix an interior edge $e\in\Gamma_{h}$ that is shared by elements $K_{1}$ and $K_{2}$ and denote by ${\bm{\Lambda}}_{K_{1}K_{2}}$ the set of degrees of freedom for the unknown $\lambda_{h}$, that lies on the edge $e$.

[TABLE]

If the edge $e$ lies on the boundary $\partial\Omega\cap\partial K$ we have

[TABLE]

Now we summarize a general HDG assembly and solve procedure:

1. (i)

Use equation (12) to assemble the local problems for $\bm{\Lambda}_{K}$. 2. (ii)

Assemble the local problems for $\bm{\Lambda}_{K}$ into a global system matrix using (13). 3. (iii)

Solve the global system matrix for $\Lambda$. 4. (iv)

Using the newly solved for $\Lambda$ in equation (13), reconstruct $\bm{U}$ and $\bm{Q}$. 5. (v)

Postprocess $U$ and $Q$ to obtain superconvergence.

Step (ii) in our framework is technically assembly free, we do not directly assemble and store the global system matrix. Instead, we exploit the unassembled local problems from equation (12) to generate a matrix vector multiplication routine. The resulting routine utilizes many small but dense matrices, instead of a large sparse matrix. See [54] for a survey of different finite element assembly strategies (matrix free or otherwise). Although the focus of our paper is not assembly, an algorithm was developed in [36] to accelerate the assembly of the HDG global system matrix.

4 Basis functions

We select a tensor product basis for the space $\mathbb{Q}^{p}(K)$. To easily facilitate high order approximations, we invoke a nodal basis with nodes that correspond to roots of Jacobi polynomials (Gauss-Legendre-Lobotto (GLL)). Various operators in the HDG discretization and multigrid method may require evaluation of the nodal basis at points that are not nodal. Some authors use a modal–nodal transformation to deal with such evaluations (see [35]). Another technique is to simply use the fast and stable barycentric interpolation (see [10]). This allows one to stay in the Lagrange basis and not have to resort to a generalized Vandermonde matrix. That is, barycentric interpolation allows for the stable evaluation of the Lagrange basis at any point in its domain. Given $N$ grid points, the setup cost is $\mathcal{O}(N^{2})$ to generate the barycentric weights, and a $\mathcal{O}(N)$ cost for each evaluation. Due do the discontinuous nature of the HDG approximation, evaluations occur element-wise, so the stability of barycentric interpolation is perhaps more important than cost of revaluation. Since barycentric interpolation allows for arbitrary evaluation, one can select quadrature points that differ from the nodal interpolation points. Below we introduce the barycentric interpolation formula used in this work.

Let $v$ be a polynomial of degree $n$, that interpolates the scattered data $\{f_{0},f_{1},\ldots,f_{n}\}$ through the points $\{x_{0},x_{1},\ldots,x_{n}\}$. The second form of the barycentric formula is

[TABLE]

where $W_{j}$ are the barycentric weights. Interpolation points $x_{i}$ arising from classical orthogonal polynomials have explicit formulas for their barycentric weights ([52], [55]). As an example, GLL and GL nodes have the following barycentric weights:

[TABLE]

where $w_{j}$ and $x_{j}$ (with appropriate superscripts) are the GL and GLL quadrature nodes and weights. On the reference element, nodal basis functions are defined using a tensor product of a 1D spectral grid. In the case of GLL nodes, $\{\xi_{i}^{(p)}\}_{i=0}^{p}$ are zeros of a particular family of Jacobi polynomials (see [35]), where $-1\leq\xi_{i}\leq 1$. We adopt the standard that the reference element in is $[-1,1]$ in 1D, $[-1,1]\times[-1,1]$ in 2D, and $[-1,1]\times[-1,1]\times[-1,1]$ in 3D. See Fig. 1 for a sample spectral grid.

In 2D the tensor product can be written as

[TABLE]

where $p$ is the polynomial degree associated with $Q^{p}$. On each element, we define the basis functions as tensor products:

[TABLE]

with $I=i+j(P+1)$ (lexicographic ordering, and $\ell_{i}$ is the Lagrange polynomial that is nodal at GLL node $i$). In one dimension, if there are $p+1$ GLL nodes, then there are $p+1$ associated basis functions. In two dimensions, we take the tensor product of the 1D GLL nodes, which results in $(p+1)^{2}$ GLL nodes, and $(p+1)^{2}$ associated basis functions. In the index $I$, if $i=j=p$, then $I=p+p(p+1)=p^{2}+2p$, and including the index contribution from $i=j=0$, we have a total of $p^{2}+2p+1=(p+1)^{2}$ basis functions in 2D as well. For brevity, let the lexicographic ordering $I=i+j(p+1)=:(i,j)$. Let the reference element be given by $\Box=[-1,1]\times[-1,1]$. Just as in the 1D case, we assume that

[TABLE]

Using a change of variables to map from the physical element to the reference element we have

[TABLE]

where $w_{i},w_{j}$ are the GLL quadrature weights in the $x$ and $y$ directions. We want to map an arbitrary element to this element, which requires a change of variables $\xi=\xi(x,y),$ $\eta=\eta(x,y)$. The Jacobian of this mapping is

[TABLE]

Hence, integrating over an arbitrary element $\Omega^{e}$ ($\cup_{e=1}^{E}\Omega^{e}$),

[TABLE]

where $\mathcal{J}^{e}(\xi_{i},\eta_{j})=\mathcal{J}^{e}(\xi,\eta)$. The above approach assumed that the interpolation points and quadrature points were the same (classic spectral element method). However, GLL quadrature rule for $p+1$ points is only exact for polynomials of degree $2p-1$. If higher order quadrature is needed, the GL quadrature rule for $p+1$ points is exact for polynomials of degree $2p+1$. Then, one can fix the GLL interpolation points, and evaluate the Lagrange basis functions at GL quadrature points if desired. To do this, we utilize the barycentric formula:

[TABLE]

Barycentric interpolation enables the stable and fast evaluation of the Lagrange polynomial basis $\ell_{j}$ anywhere in its domain.

5 HDG discretization

Before proceeding with the results of the HDG GMG method, we verify numerically that the HDG discretization provides the expected optimal $L^{2}$ convergence rates; $p+1$ for both the potential $u_{h}$ and its flux ${\bm{q}}_{h}$. Moreover, with the use of a local postprocessing filter [18], we can achieve superconvergence of the potential $u_{h}$, so that it converges in the $L^{2}$ norm with the rate $p+2$. HDG does fall under the umbrella of stabilized DG methods, so the parameter $\tau$ in equations (8) and (9) needs to be specified. The local stability parameter is piecewise constant, defined facet-by-facet. For the model problem we set $\Omega=[0,1]\times[0,1]$, $\partial\Omega_{D}=\partial\Omega$ ($\partial\Omega_{N}=\emptyset$), ${\bf K}(x,y)=\tanh{(x+y)}+1$. The domain $\Omega$ is partitioned into $N\times N$ squares. A manufactured solution is used to examine the error: $u(x,y)=x(x-1)y(y-1)\exp(-x^{2}-y^{2})$, and the corresponding forcing function $f$ is determined from $u$.

For the model Poisson problem it turns out that $\tau\equiv 1$ on every facet provides optimal convergence rates. In Table (1), it is apparent that the expected convergence rates are met. Moreover, the postprocessed potential $u_{h}^{*}$ results in a rate of $p+2$. The postprocessed flux ${\bm{q}}_{h}^{*}$ converges in a rate of $p+1$, but the errors are smaller than that of ${\bm{q}}_{h}$.

6 Multigrid algorithm

Multigrid methods work by eliminating a wide range of frequencies from the estimated error between the exact solution, and the multigrid iterate. Any multigrid method has three key components: relaxation, grid transfers, and coarse grid operators (course grid correction). In the following subsections our selection of these components is presented. The geometric multigrid method leverages information about the underlying geometry, discretization, and PDE. This results in an optimal solver ([53], [12]) that is coupled to the discretization; only $\mathcal{O}(N)$ floating point operations are required to reduce the error to discretization level accuracy. Moreover, GMG offers a uniform convergence rate. That is, the rate of convergence of GMG is not dependent on the mesh size.

In direct contrast to GMG, algebraic (black box) solvers rely entirely on the discretization matrix. The trade off with algebraic solvers (e.g. Krylov methods) is that they are decoupled from the discretization, but their convergence is very sensitive to the condition number of the discretization matrix. Refinement in $h$ or $p$ increases said condition number, and in turn the algebraic solver requires even more iterations to reach a desired tolerance. In this context, very efficient preconditioners are needed to assist the algebraic solvers. Multigrid methods are also very popular as preconditioners ([23]).

The finite element method brings with it two ways to potentially increase accuracy: $h$ refinement and $p$ refinement. Since GMG is coupled to the discretization, it can leverage $h$ refinement and/or $p$ refinement (refinement in $p$ means increasing/decreasing the polynomial order). This implies that two different classes of multigrid operators are needed. One for $h$ refinement, the other for $p$ refinement. For simplicity, we do not consider $hp$–GMG (refining $h$ and $p$ simultaneously). In $h$–GMG the polynomial order $p$ is fixed, and the hierarchy of grid spaces are determined by increasing/decreasing the mesh size. In $p$–GMG the mesh size $h$ is fixed, and the hierarchy of grid spaces are determined by their increasing/decreasing the polynomial order. The combination of these three refinement strategies offers interesting and potentially efficient ways to deal with high order and very high order approximations ([3], [22]).

We use a nested multigrid strategy. The multigrid method starts with $p$-multigrid, which fixes the mesh and the fine and coarse grids are formed by increasing/decreasing the polynomial order. Once the polynomial order reaches $p=0$, which corresponds to a face-centered finite volume scheme, we continue with the $h$-multigrid scheme as suggested in [24], [19], [tan2009iterative]. This technique exhibits $h$-independent convergence. It projects the HDG $p=0$ approximation into an appropriate space of continuous functions, so that standard $h$-multigrid for continuous Galkerin can be used for the coarser levels. A diagram of our general nested multigrid strategy is given in Fig. 2.

We use rediscretization for $p$-multigrid to form the hierarchical grids; this way we do not store the sparse linear systems for the polynomial orders $p>0$. The corresponding grid transfer operators are also implemented in a matrix free manner, as the canonical grid transfer operators are discontinuous across element interfaces, they are completely data parallel.

For $p=0$, and subsequent $h$-multigrid levels, we store the linear systems in compressed row storage format. On the coarsest level (matrix of about dimension 2000) of the $h$-multigrid scheme, we employ the MKL optimized PARDISO direct solver [48]. Other options are available for the coarse grid solve, for instance replacing the direct solve with a number of relaxation steps, or continuing the multigrid hierarchy with an algebraic multigrid method.

For problems where geometric multigrid is not suitable, or generating a meaningful geometric hierarchy is not possible, one can replace $h$-multigrid with algebraic multigrid or an appropriate coarse grid solver. In these cases, $p$-multigrid is still applicable as the mesh is fixed.

6.0.1 $p$-multigrid grid transfer

Finite element methods pair well with multilevel methods, due to the flexibility in chose of underlying spaces. Consider the nested discontinuous finite element spaces $\Omega^{H}\subset\Omega^{h}$. The mapping from a coarse space ($\Omega^{H}$) to a fine space ($\Omega^{h}$) is called prolongation, or interpolation. There are many ways to build a prolongation mapping, however, in the context of finite element methods there is a canonical choice.

The canonical prolongation operator uses the fact that $\Omega^{H}\subset\Omega^{h}$. That is, any coarse grid function can be expanded as a linear combination of fine grid functions. Thus, we take the prolongation operator to be the natural embedding from $\Omega^{H}$ into $\Omega^{h}$. As such, $P:\Omega^{H}\to\Omega^{h},$

[TABLE]

The mapping from a fine space to a coarse space is called restriction, or sub-sampling. It is not as trivial to select a restriction operator. However, since finite element methods are based on a weak formulation, a canonical restriction is defined by a $L^{2}$ projection. The restriction is given by $R:(\Omega^{h})^{\prime}\to(\Omega^{H})^{\prime}$,

[TABLE]

Notice that $(Rv_{h},u_{H})=(v_{h},R^{*}u_{H})=(v_{h},Pu_{H})$. One can infer from this that $R^{*}=P$. Selecting the restriction operator in this way is deliberate; it is done to ensure that the multigrid operator is symmetric. This allows the multigrid method to be a suitable preconditioner for symmetric iterative methods.

The grid transfer operators are defined on the reference element, and for discontinuous Galkerin methods the canonical prolongation and restriction operators are data independent. For $p$-multigrid we apply the grid transfer operators in a matrix free manner. These grid transfer operators are defined in equations (14) and (15), where it is understood that for $p\geq 0$, the coarse space is $\Omega^{H}={M}_{h}^{p}$, and the fine space is $\Omega^{h}={M}_{h}^{p+1}$.

We note that for the HDG method the grid transfer operators are defined on a lower dimension. This is due to static condensation reducing the unknowns (dimension $d$) to the mesh skeleton (dimension $d-1$). For instance, classical DG methods in 2D would have a prolongation operator $P:Q^{p}\to Q^{p+1}$ of size $(p+2)^{2}\times(p+1)^{2}$. For HDG in 2D, the prolongation operator would be $P$ of size $(p+2)\times(p+1)$, as static condensation reduces the unknowns to 1D. Nodal basis functions give rise to dense grid transfer operators, whereas modal basis functions have sparse binary grid transfer operators.

6.0.2 $h$-multigrid grid transfer

The $h$-multigrid scheme we use is taken from [24], [19], [tan2009iterative]. The main idea is to project the $p=0$ HDG solution into the space of continuous functions, then continue with a standard geometric multigrid for piecewise linear continuous finite elements. It is also possible to project the $p=1$ HDG solution to the piecewise linear finite element space. However, in [51] and [Cockburn17102013], it is found that this leads to larger iteration counts and poor convergence rates compared to the $p=0$ projection. Our numerical experiments also agree with their findings.

Let $\mathcal{E}_{k}$ denote a uniform mesh of quadrilateral elements at level $k$. The coarsest mesh is denoted by $k=0$. Geometric refinements of the coarsest mesh are then made, for a total of $J$ meshes. The mesh for the original discretization (finest mesh) is denoted by $k=J-1$. For simplicity we use uniform meshes so that a geometric hierarchy can be used. Set

[TABLE]

for $0\leq k\leq J-1$, and $M_{J}=M_{h}^{0}$. Note that the spaces have a non-nested property $M_{0}\subset M_{1}\subset\ldots\subset M_{J-1}\not\subset M_{J}.$ As such, a specialized prolongation operator has to be used [24]. It is defined as $P_{k}:M_{k-1}\to M_{k}$ ($2\leq k\leq J$),

[TABLE]

where $\Pi_{M_{J}}$ is the $L^{2}$-orthogonal projection onto $M_{J}$. The corresponding restriction operator $R_{k-1}:M_{k}\to M_{k-1}$ follows from equation (14) ($1\leq k\leq J-1$),

[TABLE]

6.0.3 Coarse grid operator

The multigrid method utilizes a hierarchy of grids. As such, there is a need for a hierarchy of discretizations. In other words, on a coarse grid we have a reduced model of our original discretization. This reduced model needs an appropriate discretization (coarse grid operator) for the coarse grid. Two of the most common coarse grid operators are subspace inheritance and subspace non-inheritance.

Subspace inheritance defines coarse grid operators by ${\bf A}^{H}={\bf R}{\bf A}^{h}{\bf P}$, where ${\bf A}^{H}$ is the coarse grid operator, ${\bf A}^{h}$ is the fine grid operators, ${\bf R}$ and ${\bf P}$ are the restriction and prolongation operators, respectively. On each coarse grid information from the previous level is being inherited. This technique has the disadvantage that continuity may be implicitly enforced for sufficiently coarse grids. Subspace non-inheritance simply uses rediscretization, on each coarse grid $\Omega^{H}$ the discretization operator is built for the dimension associated with $\Omega^{H}$.

In the $p$-multigrid setting, the coarse grid problems are formed by discretizating our problem for smaller polynomial degrees. With $p$-multigrid there is a subtle choice that is to be made for its refinement strategy. The coarse spaces for $p$-multigrid are obtained by decreasing the polynomial order $p$. How one should decrease the polynomial order is not well agreed upon. A couple of coarsening ratios have been proposed in the context of continuous $hp$-multigrid (see [43]):

• •

geometric: assumes $p$ is a power of two, $p,p/2,p/4,\ldots,2,1,2,4,\ldots$,

• •

odd: assumes $p$ is a an odd natural number,

• •

even: assumes $p$ is a an even natural number,

• •

full arithmetic: $p,p-1,p-2,\ldots,2,1,0,1,2,\ldots$.

The most flexible coarsening ratio is full arithmetic, which allows for arbitrary polynomial degrees.

To define the $p$-multigrid hierarchy we use full arithmetic coarsening. We selected full arithmetic coarsening since we are using multigrid as a solver, not as a preconditioner. In this case, there are more opportunities to eliminate high frequency modes. For high order DG methods, traditional parallelizable relaxation schemes do not properly eliminate high frequency modes [33, 34, 23, 4, 14, 13]. Full arithmetic coarsening was found to be as effective as more rapid coarsening for problems with moderate polynomial degrees [43, 22].

For $p>0$, we employ a matrix free procedure. In the $h$-multigrid setting, the coarse grid problems are formed by subspace inheritance. For low order problems, there is little benefit for matrix free operations, so we store the $p=0$ matrix and any subsequent matrices that arise from the $h$-multigrid algorithm. This amounts to the Galerkin triple product ${\bf A}_{k-1}={\bf R}_{k-1}{\bf A}_{k}{\bf P}_{k}$ (for $1\leq k\leq J-1$), where ${\bf R}_{k-1}$ and ${\bf P}_{k}$ are as defined in equations (17) and (16), respectively.

6.0.4 Relaxation

Relaxation (also called smoothing) plays an integral part in multigrid techniques. High frequencies are damped using a relaxation method, leaving low frequencies to be resolved by coarse grid correction. In some sense, multigrid can be thought of as a way to accelerate the convergence of a relaxation method; typically the relaxation methods used in multigrid have poor (or no) convergence properties. The main purpose of a relaxation method is to smooth the error, not necessarily make it small. Standard choices for relaxation are pointwise/block stationary iterative methods (Jacobi, damped Jacobi, Gauss-Seidel, Successive over relaxation, etc.), and polynomial type preconditioners.

A general overview of iterative methods can be found in [46]. For lower order discretizations, (in the context of this work, third order or less) standard relaxation methods tend to work well. However, for higher order discretizations, loss of uniform convergence is experienced. A larger number of relaxation steps, or more powerful smoothers are required to retain uniform convergence in the high order case.

In a parallel setting, polynomial smoothers where shown in [1] to perform better in terms of time-to-solution. With the advent of many-core devices, some authors have examined parallel asynchronous iterations (see [2], [5]). Sparse approximate inverses (SPAI) have exhibited promising results when used as preconditioners and smoothers in multigrid (see [9], [15]). Moreover, the construction of SPAI operators are inherently parallel, and only require a matrix vector product to apply.

Relaxation for high order discontinuous Galkerin methods typically require specialized treatment. This was found to be the case for HDG in our experiments, as well as in [tan2009iterative] and [58]. In particular, block-type methods work better than pointwise or polynomial smoothers. In the context of our work, we have found that SPAI and additive Schwarz domain decomposition methods work very well as a smoother. In Section 7 we present multigrid convergence results for an SPAI preconditioner. The additive Schwarz smoother has similar convergence properties depending on the overlap size. These more expensive smoothers are used in our work as they able to retain robust multigrid convergence. Section 7 discusses this matter further.

6.0.5 Sparse approximate inverse relaxation

Sparse approximate inverses have exhibited promising results when used as preconditioners and smoothers in multigrid (see [9], [15]). Moreover, the construction of SPAI operators are inherently parallel, and only require a matrix vector product to apply. SPAI seeks an approximation ${\bm{M}}$ of ${\bm{A}}^{-1}$ such that

[TABLE]

where $\mathcal{S}$ is a set of sparse matrices. By using the Frobenius norm, a great deal of concurrency is exposed:

[TABLE]

Each of the $n$ least squares problems are independent of one another. Storage and computational savings can be obtained by leveraging the sparsity of $\bm{A}$ and $\bm{M}$. In addition, since the HDG method gives rise to a symmetric operator, the factorized sparse approximate inverse (FSAI) can be utilized to obtain further improvements. The use of SPAI as a smoother will require its explicit storage, but the application is that of a matrix vector multiply:

[TABLE]

6.0.6 Additive Schwarz domain decomposition relaxation

Increasing the polynomial order in finite element methods causes the condition number of the discretization operator to grow. This stretches the spectrum of the discretization operator, and traditional smoothers can no longer damp a sufficiently wide range of high frequencies. To combat this, one can increase the number of smoothing steps, or seek a more powerful smoother. In some cases, domain decomposition methods can provide a very powerful smoother for multigrid techniques. In [41], [49], and [50], variations of the Schwarz domain decomposition method is used as a multigrid smoother for high order finite element discretizations.

The standard additive Schwarz method takes the form

[TABLE]

where ${\bm{R}}_{e}$ is a binary restriction operator that transfers global datum to local datum. The operator ${\bm{A}}_{e}$ is the discretization operator restricted on the facet $e$. Notice that this is a single level Schwarz method, and as such, it is not as effective as their two level counter parts. Weighting the single level Schwarz method has been demonstrated to produce effective smoothers (see [41],[50], [40]). The weighting appears to be successful, but is also somewhat experimental. The inverse of a diagonal counting matrix is used in [41], and a type of bump function is used in [50] (both of these weightings are used without a convergence theory). For a weighting matrix $\bm{W}$, the additive Schwarz smoother is applied as follows:

[TABLE]

On each level of the $p$-multigrid hierarchy we store smoother $\bf M$ as dense blocks, and apply them in a unassembled manner. The size of these blocks depends on the polynomial order and overlap size.

6.0.7 Standard multigrid V-cycle

Our nested multigrid V-cycle is presented in Algorithms 1 and 2. This is often referred to as the $V$-cycle. We use the following notation: $S^{\nu}(\cdot)$ denotes $\nu$ steps of a relaxation method. In the $p$-multigrid cycle the Smoother $S^{\nu}(\cdot)$ is either SPAI or the additive Schwarz relaxation. In Section 7 we use SPAI, and in Section 8 we use the additive Schwarz relaxation. For $h$-multigrid we use the Chebyshev smoother [46] as the problem is low order.

7 Multigrid convergence

In [3], uniform convergence in $h$ and $p$ was established for a discontinuous Galerkin $hp$-multigrid algorithm. However, in order to guarantee $p$ (or $h$) independent convergence, the number of relaxation steps per level grows quadratically in $p$. Numerical experiments suggest that for our HDG $p$-multigrid, the relaxation steps per level need to grow linearly (in $p$) to guarantee $p$ independent convergence. In [19], $h$ independent convergence was proven for the HDG $h$-multigrid algorithm we use.

Fig. 3a and Fig. 3b display the results for $2\leq p\leq 8$ and fine mesh with $(2^{5})^{2}$ elements. We employ a FSAI smoother on each level, with $\nu_{1}=\nu_{2}=2$ pre and post smoothing steps. The FSAI smoother is constructed so that it results in an approximate inverse with a operator complexity of unity (nnz[$\bm{M}$]$=$nnz[$\bm{A}$]). Subspace non-inheritance is used to generate coarse grid operators. For the $p$-GMG phase of the GMG method, we use full arithmetic coarsening. To measure the convergence rate, we keep track of the two norm of the fine grid residual between successive iterations:

[TABLE]

We can see that the results are quite good, even for a modest FSAI smoother. All of the convergence rates are under 0.22, and tend to cluster in the range 0.13-0.15 for even polynomial degrees. This observation perhaps indicates that even-$p$ coarsening or geometric-$p$ will not only be more efficient (more rapid coarsening), but also more accurate.

Fig. 3 seems to suggest that our method is not exhibiting a $p$ independent convergence rate. However, we recall that according to [3], uniform convergence in $p$ guaranteed only if the number of relaxation steps per level grows quadratically in $p$ (we set $\nu_{1}=\nu_{2}=2$). For multigrid methods there is a subtle observation to be made: uniform convergence is an excellent property; but if the uniform convergence rate is near one, this is not an efficient solver. Moreover, if the number of smoothing steps is sufficiently large, a great deal of efficiency is lost. For reasonable smoothing steps, in [14] and [3] the authors demonstrate convergence rates of 0.85 (or worse). In order to drop the convergence rates down to $0.5$, upwards of 20 smoothing steps every level needs to be taken.

What is needed for efficiency is as follows: uniform convergence, a fast convergence rate (problem dependent), and a small number of uniform smoothing steps. Satisfying these three properties in practice for high order discontinuous Galkerin methods is not trivial. Some potential avenues to explore are stronger (but more expensive) smoothers, different coarse grid transfer operators, different multigrid schedules, as well as different grid hierarchies. This direction of inquiry is beyond the scope of our paper. In Fig. 4 we do show that a more expensive smoother essentially yields uniform convergence with a fast convergence rate and a fixed number of smoothing steps.

For the next set of experiments, we allow the FSAI operator complexity to grow ((nnz[$\bm{M}$]$\,\approx\,2.5\,$nnz[$\bm{A}$]). In Figs. 4a, 4b, 4c, and 4d, display the results of this change. For the $V$-cycle, the residual is reduced to machine precision after only 5 to 7 iterations (Fig. 4a). The aggressive smoothing also yields stellar convergence rates, below 0.015, as can be ascertained from Fig. 4b.

The performance of the V-cycle is better than ideal, so we can easily extend our HDG GMG method to leverage the full multigrid cycle (FMG). FMG is well known to be the optimal multigrid schedule for linear problems; with a single FMG iteration, the residual is reduced to discretization level error. Moreover, it only requires $\mathcal{O}(N)$ floating point operations to achieve this accuracy ([11], [12], [53]). Indeed, Fig. 4c numerically verifies that a single FMG iteration is enough to reach discretization level error. There is something that particularly interesting about Fig. 4c - the FMG iteration performs better for higher orders $6\leq p$. This second experiment of course comes at a price: the FSAI smoothers on each level have an operator complexity in the rage of 2.65 to 2.85. Such a trade off may be valuable for problems with highly varying or discontinuous coefficients. Also, the FSAI smoother can be easily tuned to control how aggressively its operator complexity grows.

8 Performance model

The storage and assembly of global matrices in finite element methods can be exceedingly prohibitive, especially at higher orders. By leveraging matrix free algorithms, one can save on memory, and, convert a memory bound problem (sparse matrix vector multiplication) into a compute bound problem. The authors in [45] found that to improve sparse matrix vector multiplication for HDG methods, specialized storage formats were needed. Since the HDG method can be reduced to a problem on the trace space, this allows for a assembly (matrix free or otherwise) of the discretization operator in a facet-by-facet manner, instead of a element-by-element manner. The importance of this is that the element-by-element approach requires a synchronization (barrier, atomic, etc.) in order to avoid race conditions. A graph coloring algorithm is typically used in this situation, but only allows for a group of colors to be utilized at any given time. The facet-by-facet approach does not depend on the vertex degree, but only on the element type (how many facets on a given element). Fig. 5 displays some example DOFs and connectivity for HDG, DG, and continuous Galerkin methods.

The HDG method is able to exploit facet-by-facet connectivity; given an interior facet, only the two elements that share said facet will contribute to its DOFs. An element-to-facet mapping will allow one to gather and scatter the DOFs on a facet. Pseudocode for our matrix vector multiply routine is given in Algorithm (3).

8.0.1 Knights Landing (KNL) many-core coprocessor

In this work we use the second generation Xeon Phi designed by Intel®. It is a many-core processor, and the 7000 series has anywhere from 64 to 72 cores, with 4 threads per core. The KNL is an example of a high throughput low memory device. The design of the KNL is similar to other many-core devices and accelerators: a large number of cores with lower clock speeds are packed into the unit, enabling a large vector width, as well as having access to a user-manageable fast memory hierarchy. One interesting feature of the KNL is that it can be programmed using traditional parallel paradigms like OpenMP, MPI, and pthreads. In addition, it operates as a native processor. Further details regarding the KNL can be obtained in [32]. For specifications of the KNL used in this work, see Table 2.

8.0.2 Local matrix generation

In order to make use of the assembly free matrix vector multiply, one needs to generate the associated local matrices (in the case of HDG, see (i), (iv), and (v) in Section 3). To generate the required local matrices, for simplicity, we use a one core (thread) per element strategy. However, as one increases the polynomial order (beyond $p=5$), this strategy loses performance. Further improvements may potentially be obtained by utilizing nested parallelism, or linear algebra libraries dedicated for small or medium sized matrices ([26], [36]).

For all computational experiments we set the KNL in quadrant mode. Figs. 7 and 7 display the wall-clock and speed up as the polynomial order and number of threads is varied. In Fig. 7, one can see that as we increase the number of threads beyond 32 or 64, diminishing returns are very noticeable. To generate the results of Fig. 7, we fix the number of threads to 64, and we vary the polynomial order. The comparison with a serial implementation clearly shows that the additional parallelism the KNL offers is beneficial; speed ups ranging from 2X to 32X are attained.

8.1 Roofline analysis

The KNL (7210 processor number) used in this analysis runs at 2.1 GHz (double precision), a double precision processing power of of 2,199 GFLOP/s, and the STREAM memory bandwidth benchmark (triad, see [42]) reports a bandwidth of 300 GB/s. These two metrics provide the performance boundaries for arithmetic throughput and memory bandwidth limits, respectively. For our numerical experiments, the clustering mode is set as quadrant, and cache mode is set to flat.

8.1.1 Sparse Matrix Vector Multiply (SPMV)

In general, matrix vector multiplication is severely bandwidth bound. For a dense matrix of dimension $N$, we can expect $(2N^{2}-N)$ FLOPs and $8(2N^{2}+N)$ MEMOPs (bytes). Thus, for large $N$, we can expect an arithmetic intensity of only $\lim_{N\to\infty}(2N^{2}-N)/(8(2N^{2}+N))=0.25$. For finite element problems the underlying discretization matrix is typically sparse, which pushes SPMV further into the bandwidth bound region on the roof line chart. Matrix free operation can remedy this situation by improving the arithmetic intensity. That is, matrix free application can shift our bandwidth bound SPMV to a compute bound (or less bandwidth bound) problem.

In Fig. 9a, we include a roofline analysis of our matrix vector multiply routine; in the context of its performance on the Laplacian operator. The roofline analysis ([57]) allows us to identify bottlenecks, confirm algorithm limitations, as well as gives us insight on what we should focus on in terms of optimization. Two different techniques are studied. Sparse matrix storage (CSR format, see Fig. 9b) has the lowest performance arithmetic intensity, with a theoretical limit of 0.25. The multithreaded Intel MKL library is used for this approach ([31]). For high orders the HDG method has hundreds of nonzeros per row. Slightly better but similar performance was obtained in [45] by using a specialized block sparse matrix storage format on GPUs. The matrix free technique shifts the arithmetic intensity favorably (see 9a); this behavior is typical in spectral methods and spectral elements, due to near constant FLOP and memory requirements per degree of freedom ([54], [16], [39]).

Fig. 8 displays the bandwidth measurements. From Fig. 8a, for lower order polynomials ($p<5$) on a fixed coarse mesh, we see that a maximum of 50% of the peak bandwidth is achieved. The overhead costs of forking threads is dominating computations for $p<5$ on coarser meshes. As $p$ increases beyond five, the cost of forking threads is no longer dominant. Our algorithm reaches 80% of the peak bandwidth reported by the STREAM benchmark.

Fig. 8b considers fixing the polynomial degree, and increasing the number of elements in the mesh. Here it is evident that $h$-refinement improves the bandwidth performance for lower polynomial degrees. To observe further improvements, a new SPMV routine could be designed specifically for lower orders. For instance, in [38], the authors group together the contributions from multiple elements into a thread block in their SPMV-GPU kernel. Their technique is shown to be very high performing.

8.1.2 Projection From Volume To Surface

The generation of the local solvers from equation (12) is completely data parallel and needs no synchronizations or thread communication. It requires a dense linear solve, two matrix vector multiplications and a single SAXPY. Roughly the complexity model for this projection operator in terms of FLOPs can be estimated by

[TABLE]

where $N=4(p+1)^{2}$ and $M=4(k+1)$. Table (3) collects the results for polynomial order $0,1,2,3,$ and $4$. As we increase the polynomial order, performance increases due to the compute bound nature of the local HDG solvers; as the dense linear solve dominates with $\mathcal{O}(n^{3})$ FLOPs for $\mathcal{O}(n^{2})$ data. We again note that the focus of this work is not on the generation of this local matrices, but we can still obtain reasonable performance utilizing a straightforward implementation.

After solving the trace space system given by equation (13), if the volume solutions $u$ and $\bm{q}$ are desired, one can invoke equation (12) to reconstruct them. The cost of this procedure depends on if one discards the local matrices. In this case, one has to recompute and the cost is the same as the projection from the surface to volume. If one instead keeps the local matrices, all that is needed is two matrix vector multiplies and a single SAXPY. Ultimately this comes down to a preference of convenience over memory concerns.

8.1.3 Cost analysis matrix-free vs matrix-stored

We consider the cost of our matrix free approach compared to a matrix based approach. The previous sections discuss the memory considerations, but there are other measures one can use. Most importantly, operation counts and time-to-solution. The operation count (FLOPs) for sparse matrix vector multiplication is $\mathcal{O}(nnz)$, where $nnz$ is the number of nonzero entries in the sparse matrix in question. For the statically condensed HDG method on a structured Cartesian mesh, $nnz$ for the stiffness matrix is $C_{1}=\mathcal{O}(Nd(p+1)^{d-1}\times(4d-1)(p+1)^{d-1})$ [47], [30] ($d$ is the dimension, $N$ is the total number of elements in the mesh). The operation count for our proposed matrix free variant (Algorithm 3) is $C_{2}=\mathcal{O}(2|\Gamma_{h}|(2\cdot 4(p+1)-1))$. In 2D, as $N\to\infty$, the ratio tends to $C_{2}/C_{1}\to(16p+14)/(7(p+1)^{2})$. Which indicates that for higher orders, the matrix free variant is cheaper in terms of operations than the sparse matrix format.

A second important measure is time-to-solution. Both the matrix free and matrix based algorithm require the local solvers. In the matrix based approach, a non negligible amount of time is spent on the matrix assembly, especially for higher orders. We compare the wall-clock matrix based solve to the wall-clock of the (parallel) assembly. See Fig. 10a. For $p=2$, the assembly time is already more than 40% of the total solve time. Increasing the polynomial order to $p>2$, the assembly time ranges from around 45% to 56% of the total solve time. Comparing the matrix free solve timings, to the matrix based solve and assembly time, Fig. 10b shows for $p>2$ we obtain speedups of 1.5 to 2.3 in favor of the matrix free method.

We acknowledge that the benefits of matrix free methods will vanish if efficient preconditioners are not available. In such cases, in order to have more options for robust preconditioners, the matrix will have to be stored.

9 Conclusions

We presented a highly efficient multigrid technique for solving a second order elliptic PDE with the HDG discretization. Both the multigrid technique and the HDG discretization are amenable to high throughput low memory environments like that provided by the KNL architecture; this was demonstrated by using thorough profiling and utilizing a roofline analysis ([57]). The HDG method has much of its computation localized due to the discontinuous nature of the solution technique. Moreover, the HDG method brings static condensation to DG methods, which significantly reduces the number of nonzeros in the discretization operator. This in turn means that less work is required for a linear solver to obtain solutions. Since the HDG method converges with optimal orders for all of the variables it approximates, a local element by element postprocessing is available.

Any iterative solver will require sparse matrix vector multiplication (matrix free or otherwise); including multigrid. Two different approaches were examined for sparse matrix vector multiplication: matrix based and matrix free. For high order HDG, the matrix free technique better utilizes the resources available on the KNL, because of increased arithmetic intensity. In the high order regime, the matrix based technique requires a sparse matrix storage which increases time to solution. This is mainly because of additional assembly time, having a low arithmetic intensity, and erratic memory access patterns for sparse matrix vector multiplication.

Our algorithm is able to attain 80% of peak bandwidth performance for higher order polynomials. This is possible due to the data locality inherent in the HDG method. Very good performance is realized for high order schemes, due to good arithmetic intensity, which declines as the order is reduced. The performance is lower for polynomial orders $p<6$. With an approach similar to the work done in [38], where multiple cells are processed by a thread block, one would be able to increase performance for lower order polynomials. We observed speedups when compared to a multicore CPU for the HDG methods components, namely, volume to surface mapping, surface to volume mapping, matrix free matrix vector multiplication, and local matrix generation. The ratio of peak flop rates on the two target architectures was roughly 100X, and peak bandwidths was 5X, so this figure fits with our model of the computation. This is possible due to the data locality inherent in the HDG method. Very good performance is realized for high order schemes, due to good arithmetic intensity, which declines as the order is reduced.

An attractive feature of the KNL is that it can be programmed with traditional parallel paradigms like OpenMP, MPI, and pthreads. Thus, one can harness the massive fine-grain parallelism that the KNL offers by utilizing these traditional parallel paradigms with significantly limited intrusion.

Acknowledgments

Fabien acknowledges the support from the Ken Kennedy Institute and the Ken Kennedy–Cray graduate Fellowship. Fabien also acknowledges the support from the Richard Tapia Center for Excellence & Equity, and the Rice Graduate Education for Minorities program. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575. The KNL used in this work was provided to the RELACS research group by Richard T. Mills through the Intel® Parallel Computing Center at Rice University.