1 – DG Methods for Diffusion Problems

Per-Olof Persson

Introduction to the DG Method in 1-D

The Finite Difference Method (FDM)

  • Consider linear convection: \(u_t+u_x=0\) for \(x\in[0,1]\), \(u(0)=u(1)\)

  • Approximate \(u_x\) point-wise using difference formulas:

\[\begin{aligned} \frac{d}{dx} u(x_n) \approx \frac{u_{n+1}-u_{n-1}}{2\Delta x} \end{aligned}\]

  • or high-order:


\[\begin{aligned} \frac{d}{dx} u(x_n) \approx \frac{u_{n+2}-8u_{n+1}+8u_{n-1}-u_{n-1}}{12\Delta x} \end{aligned}\]

  • or one-sided (e.g. for stability, “upwinding”):

\[\begin{aligned} \frac{d}{dx} u(x_n) \approx \frac{3u_n-16u_{n-1}+36u_{n-2}-48u_{n-3}+25u_{n-4}}{25\Delta x} \end{aligned}\]

  • Simple, efficient, flexible

  • Needs structured neighborhood of nodes – hard to generalize to unstructured grids in 2-D and 3-D

The Finite Element Method (FEM)

  • Discretize domain into elements (intervals)

  • Seek approximate solution in space of piecewise polynomials \(\hat{X}\)

  • Impose equation weakly: Seek \(\hat{u}\in \hat{X}\) such that for all \(v\in \hat{X}\):

\[\begin{aligned} &\int_0^1 (\hat{u}_t+\hat{u}_x) v\, dx \\ &= \int_0^1 \hat{u}_t v\,dx + \int_0^1 \hat{u}_x v\, dx \\ &= \int_0^1 \hat{u}_t v\,dx - \int_0^1 \hat{u} v_x\, dx = 0 \end{aligned}\]


  • Leads to semi-discrete system \(M \boldsymbol{u}_t + K\boldsymbol{u} = 0\), with element-wise local \(M,K\) matrices

  • \(M^{-1}\) dense \(\Longrightarrow\) Explicit methods for \(\boldsymbol{u}_t = -M^{-1}K\boldsymbol{u}\) not practical

  • Also, unclear how to stabilize by upwinding (but other techniques exist, such as Streamline Upwind Petrov-Galerkin)

The Discontinuous Galerkin Method

  • Do not enforce continuity – allow “jumps” between elements



  • Galerkin formulation for single element \(\kappa=[0,h]\): For all \(v\in P^p(\kappa)\),

\[\begin{aligned} \int_0^h (\hat{u}_t+\hat{u}_x)v\, dx &= \int_0^h \hat{u}_t v\, dx + \int_0^h \hat{u}_x v\, dx \\ &= \int_0^h \hat{u}_t v\, dx - \int_0^h \hat{u}v_x\,dx +\mathcal{U}(u^+,u_p) v(h) - \mathcal{U}(u_0,u^-)v(0) \end{aligned}\]

  • Numerical flux function \(\mathcal{U}(u_R,u_L)\) allows for stabilization by high-order upwinding, e.g. \(\mathcal{U}(u_R,u_L) = u_L\)

The Discontinuous Galerkin Method

  • The DG formulation leads to linear system of equations: \[\begin{aligned} M \boldsymbol{u}_t+K\boldsymbol{u} +\begin{pmatrix} -u^- & 0 & \cdots & 0 & u_p \end{pmatrix}^T = 0 \end{aligned}\]

  • For example, with \(p=2\): \[\begin{aligned} \boldsymbol{u}_t &= -M^{-1} K \boldsymbol{u} - M^{-1} \begin{pmatrix} -u^- & 0 & u_2 \end{pmatrix}^T \\ &= \frac{1}{h} \begin{pmatrix} -6 & -4 & 1 \\ 2.5 & 0 & -1 \\ -4 & 4 & -3 \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \end{pmatrix} + \frac{1}{h} \begin{pmatrix} 9 \\ -1.5 \\ 3 \end{pmatrix} u^- \end{aligned}\]

  • Element-wise local FD-type stencil

  • Stabilized, “upwinded” through \(u^-\)

  • Extends naturally to other PDEs, N-D, unstructured meshes


The DG Scheme – Details of Discretization

  • Consider the 1-D conservation law \[\begin{aligned} \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0 \end{aligned}\]

  • Seek a solution in space of piecewise polynomial functions \(X_h\)

  • Nodal representation with values \(u_i^k\) for local node \(i\) in element \(k\): \[\begin{aligned} u_h(x) = \sum_{k=1}^n \sum_{i=0}^p u_i^k \phi_i^k(x) \end{aligned}\]

  • Example, piecewise linear functions (\(p=1\)):

image image

The DG Scheme – Details of Discretization

  • Galerkin formulation: Find \(u_h\in X_h\) such that \[\begin{aligned} \int_0^1 \frac{\partial u_h}{\partial t}v\, dx + \int_0^1 \frac{\partial f(u_h)}{\partial x}v\,dx = 0 \end{aligned}\]

  • Set \(v=\phi_i^k\) and integrate by parts \[\begin{aligned} \int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t} \phi_i^k\,dx + \left[ f(u_h(x))\phi_i^k(x) \right]_{x_{k-1}}^{x_k} -\int_{x_{k-1}}^{x_k} f(u_h)\frac{d\phi_i^k}{dx}\,dx = 0 \end{aligned}\]

  • Use a numerical flux function \(F(u_R,u_L)\) at the discontinuities \[\begin{aligned} &\int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t} \phi_i^k\,dx + F(u_0^{k+1},u_p^k)\phi_i^k(x_k) - F(u_0^k,u_p^{k-1})\phi_i^k(x_k) \\ &-\int_{x_{k-1}}^{x_k} f(u_h)\frac{d\phi_i^k}{dx}\,dx = 0 \end{aligned}\]

The DG Scheme – Details of Discretization

  • Example: \(f(u)=u\), \(F(u_R,u_L)=u_L\) \[\begin{aligned} &\hspace{-1.2cm}\int_{x_{k-1}}^{x_k} \frac{\partial}{\partial t} \left( \sum_{k=1}^n \sum_{j=0}^p u_j^k \phi_j^k(x) \right) \phi_i^k\,dx - \int_{x_{k-1}}^{x_k} \left( \sum_{k=1}^n \sum_{j=0}^p u_j^k \phi_j^k(x) \right) \frac{d\phi_i^k}{dx}\,dx \\ & +u_p^k\phi_i^k(x_k)-u_p^{k-1}\phi_i^k(x_{k-1}) = 0 \end{aligned}\]

  • Rearrange to obtain a linear system of equations \[\begin{aligned} M^k \dot{u}^k-C^k u^k + \begin{bmatrix} -u_p^{k-1} & 0 & \cdots & 0 & u_0^k \end{bmatrix}^T = 0 \end{aligned}\] for element \(k\), with elementary matrices \[\begin{aligned} M_{ij}^k = \int_{x_{k-1}}^{x_k} \phi_i^k\phi_j^k\,dx\text{ and } C_{ij}^k = \int_{x_{k-1}}^{x_k} \frac{d\phi_i^k}{dx}\phi_j^k\,dx \end{aligned}\]

Calculating Elementary Matrices

  • Consider an element of degree \(p\), width \(h\), and a nodal basis at the points \(s_i=h_i/p\), \(i=0,\ldots,p\)

    • For \(p\) high \((>4)\), use Gauss-Lobatto points instead

  • Write basis functions in monomial form \(\phi_i(s) = \sum_{j=0}^p c_i^j s^j\)

    • For \(p\) high \((>4)\), use an orthogonal basis instead

  • Nodal basis functions are defined by \[\begin{aligned} \phi_i(s_k) = \delta_{ij} = \begin{cases} 1 & i=j \\ 0 & i\ne j \end{cases} \end{aligned}\]

  • Produces a linear system of equations


Calculating Elementary Matrices

  • The linear system of equations has the form \[\begin{aligned} \begin{pmatrix} 1 & s_0 & s_0^2 & \cdots & s_0^p \\ 1 & s_1 & s_1^2 & \cdots & s_1^p \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & s_p & s_p^2 & \cdots & s_p^p \\ \end{pmatrix} \begin{pmatrix} c_0^0 & c_1^0 & \cdots & c_p^0 \\ c_0^1 & c_1^1 & \cdots & c_p^1 \\ \vdots & \vdots & \ddots & \vdots \\ c_0^p & c_1^p & \cdots & c_p^p \\ \end{pmatrix} = \begin{pmatrix} 1 & & & \\ & 1 & & \\ & & \ddots & \\ & & & 1 \end{pmatrix} \end{aligned}\] or \(VC=I\), which gives the coefficient matrix \(C=V^{-1}\)

  • Use Gaussian quadrature or explicit polynomial integration to compute the elementary matrices \[\begin{aligned} M_{ij}&=\int_0^h \phi_i(s)\phi_j(s)\, ds \\ C_{ij}&=\int_0^h \phi_i'(s)\phi_j(s)\, ds \end{aligned}\]

DG formulations for 1-D Poisson

DG for Elliptic Problems – Historical Overview

  • Enforcing Dirichlet conditions by penalties

    • Lions (1968), Babuška (1973) – Penalty term

    • Nitsche (1971) – additional terms in bilinear form for consistency

  • Interior Penalty (IP) methods

    • Babuška and Zlámal (1973) – enforce \(C^1\)-continuity by penalties

    • Wheeler (1978), Arnold (1979) – Nitsche’s method for spaces of discontinuous piecewise polynomials

  • DG methods

    • Bassi and Rebay (1997) – apply RKDG to unknown and its gradient

    • Cockburn and Shu (1998) – generalized the ideas, the LDG method

  • Unification

    • Arnold, Brezzi, Cockburn, Marini (2000,2002) – showed that most methods fit in a unified framework by choosing appropriate numerical fluxes

Second-order Equations

  • Consider the 1-D Poisson equation \[\begin{aligned} -\frac{d ^2 u}{d x ^2} = f(x)\quad\text{ in } [0,1] \end{aligned}\] with homogeneous Dirichlet conditions \(u(0)=u(1)=0\)

  • Standard Continuous Galerkin FEM would consider the space \(X_{h,0}\) of continuous piecewise polynomials satisfying the Dirichlet conditions, and solve for \(u_h \in X_{h,0}\) s.t. \[\begin{aligned} \int_0^1 -\frac{d^2u_h}{dx^2} v \,dx = \int_0^1 \frac{du_h}{dx} \frac{dv}{dx}\,dx - \left[\frac{du_h}{dx} v\right]_0^1 = \int_0^1 \frac{du_h}{dx} \frac{dv}{dx}\,dx = \int_0^1 f\,v\,dx \end{aligned}\] for all \(v\in X_{h,0}\)

  • With discontinuous functions, appropriate numerical fluxes must be chosen at all element boundaries

The 1-D Poisson Equation

  • To define a DG discretization, first split into first order system: \[\begin{aligned} -\sigma'=f(x),\quad u'=\sigma \end{aligned}\]

  • Multiply by test functions \(v,\tau\), integrate over an element, and integrate by parts to obtain the weak form \[\begin{aligned} \int_{x_k}^{x_{k+1}} f(x)\,v\,dx &= \int_{x_k}^{x_{k+1}} -\sigma' v\,dx = \int_{x_k}^{x_{k+1}} \sigma v'\,dx -\left[\hat{\sigma} v\right]_0^1 \\ \int_{x_k}^{x_{k+1}} \sigma\,\tau\,dx &= \int_{x_k}^{x_{k+1}} u'\,\tau\,dx = -\int_{x_k}^{x_{k+1}} u\,\tau'\,dx + \left[ \hat{u} \tau \right]_0^1 \end{aligned}\]

  • Galerkin formulation: Find \(u_h,\sigma_h\in X_h\) s.t. for all elements \(k\) \[\begin{aligned} \int_{x_k}^{x_{k+1}} \sigma_h v'\,dx &= \int_{x_k}^{x_{k+1}} f(x)\,v\,dx + \left[\hat{\sigma}(u_h,\sigma_h) v\right]_0^1, && \forall v\in X_h \\ \int_{x_k}^{x_{k+1}} \sigma_h\,\tau\,dx &= -\int_{x_k}^{x_{k+1}} u_h\,\tau'\,dx + \left[ \hat{u}(u_h) \tau \right]_0^1, && \forall \tau\in X_h \end{aligned}\]

  • Remains only to define the numerical fluxes \(\hat{u}(u_h),\hat{\sigma}(u_h,\sigma_h)\)

The BR1 Fluxes

  • The BR1 fluxes: \[\begin{aligned} \hat{u}=\{ u_h \},\qquad \hat{\sigma} = \{ \sigma_h \} \end{aligned}\] where \(\{\cdot\}\) is the averaging operator


  • For example, with notation according to the figure: \[\begin{aligned} \hat{u}(0) &= (u^-+u_0)/2\text{\ \ and\ \ }\hat{u}(h) = (u_3+u^+)/2 \\ \hat{\sigma}(0) &= (\sigma^-+\sigma_0)/2\text{\ \ and\ \ }\hat{\sigma}(h) = (\sigma_3+\sigma^+)/2 \end{aligned}\]

  • Simple, intuitive (no preference to direction in equation)

  • However, unstable and non-compact stencil

Interior Penalty (IP)

  • In the interior penalty method, we set \[\begin{aligned} \hat{u}&=\{ u_h \} \\ \hat{\sigma} &= \{ \nabla u_h \} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right] \end{aligned}\]

  • for some \(C_{11}>0\), where \(\{\cdot\}\) is the averaging operator and \(\left[\left[\right.\right.\cdot \left.\left.\right]\right]\) is the jump operator


  • For example, with notation according to the figure: \[\begin{aligned} \hat{u}(0) &= (u^-+u_0)/2\text{\ \ and\ \ }\hat{u}(h) = (u_3+u^+)/2 \\ \hat{\sigma}(0) &= (\left.u_h'\right|_{x=0^-} + \left.u_h'\right|_{x=0^+})/2 + C_{11} (u^- - u_0) \\ \hat{\sigma}(h) &= (\left.u_h'\right|_{x=h^-} + \left.u_h'\right|_{x=h^+})/2 + C_{11} (u_3-u^+) \end{aligned}\]

  • Convergent with optimal order of accuracy

  • However, \(C_{11}\) is problem dependent, introduces stiffness

The Local Discontinuous Galerkin (LDG) Method

  • In the LDG method, we set \[\begin{aligned} \hat{u}&=\{ u_h \} + C_{12} \left[\left[\right.\right.u_h \left.\left.\right]\right]\\ \hat{\sigma} &= \{ \sigma_h \} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]-C_{12} \left[\left[\right.\right.\sigma_h \left.\left.\right]\right] \end{aligned}\]


  • For the special cases \(C_{11}=0\) (minimal dissipation LDG) and \(C_{12}=1/2\) we get a simple upwind/downwind structure

  • For example, with notation according to the figure: \[\begin{aligned} \hat{u}(0) &= u_0 \text{\ \ and\ \ }\hat{u}(h) = u^+ \\ \hat{\sigma}(0) &= \sigma^- \text{\ \ and\ \ }\hat{\sigma}(h) = \sigma_3 \end{aligned}\]

  • Simple and general

  • Convergent with optimal order of accuracy

  • However, in general a non-compact stencil in higher dimensions

Higher Space Dimensions – Unified Framework

Higher Space Dimensions

  • From Arnold, Brezzi, Cockburn, Marini (2002)

  • Model problem: \[\begin{aligned} - \Delta u = f \quad \text{in }\Omega, \qquad u=0\quad\text{on }\partial\Omega \end{aligned}\]

  • Rewrite as first-order system \[\begin{aligned} \sigma=\nabla u, \quad -\nabla\cdot\sigma = f\quad\text{in }\Omega, \quad u=0\quad\text{on }\partial\Omega \end{aligned}\]

  • Multiply by test functions \(\tau,v\), integrate over element \(K\), integrate by parts \(\Rightarrow\) weak formulation: \[\begin{aligned} \int_K \sigma\cdot\tau\,dx &= -\int_K u \nabla\cdot\tau\,dx + \int_{\partial K} u\, n_K\cdot\tau\,ds \\ \int_K \sigma\cdot\nabla v\,dx &= \int_K f\,v\,dx + \int_{\partial K} \sigma\cdot n_K\,v\,ds \end{aligned}\]

Higher Space Dimensions

  • Introduce finite element spaces for triangulation \(\mathcal{T}_h=\{K\}\): \[\begin{aligned} V_h & := \{ v\in L^2(\Omega) : v|_K \in \mathcal{P}_p(K)\quad \forall K\in\mathcal{T}_h\ \} \\ \Sigma_h & := \{ \tau \in [L^2(\Omega)]^2 : \tau|_K \in [\mathcal{P}_p(K)]^2\quad\forall K\in\mathcal{T}_h\ \} \end{aligned}\]

  • The flux formulation: Find \(u_h\in V_h\) and \(\sigma_h \in \Sigma_h\) s.t. \[\begin{aligned} \int_K \sigma_h\cdot\tau\,dx &= -\int_K u_h \nabla\cdot\tau\,dx + \int_{\partial K} \hat{u}_K\,n_K\cdot\tau\,ds, && \forall\tau\in[\mathcal{P}_p(K)]^2 \\ \int_K \sigma_h\cdot\nabla v\,dx &= \int_K f\,v\,dx + \int_{\partial K} \hat{\sigma}_K\cdot n_K \,v\,ds, && \forall v\in \mathcal{P}_p(K) \end{aligned}\] for all elements \(K\in\mathcal{T}_h\)

  • Need to define the numerical fluxes \(\hat{u}_K\) and \(\hat{\sigma}_K\)

Higher Space Dimensions

  • Denote the union of the element edges \(\Gamma\), the interior edges \(\Gamma^0 := \Gamma \backslash \partial \Omega\), and the trace space \(T(\Gamma) := \prod_{K\in\mathcal{T}_{h}}\ L^2(\partial K)\)

  • For an interior edge \(e\), with unit normal vectors \(n_1,n_2\) define the jump and average of \(q\in T(\Gamma)\) by \[\begin{aligned} \{q\} = \frac12 (q_1+q_2) ,\quad \left[\left[\right.\right.q \left.\left.\right]\right]= q_1n_1+q_2n_2 \end{aligned}\] and for \(\sigma\in[T(\Gamma)]^2\) by \[\begin{aligned} \{\sigma\} = \frac12 (q_1+q_2) ,\quad \left[\left[\right.\right.\sigma \left.\left.\right]\right]= \sigma_1\cdot n_1+\sigma_2\cdot n_2 \end{aligned}\]

  • For boundary edges, set \[\begin{aligned} \left[\left[\right.\right.q \left.\left.\right]\right]= qn,\quad \{\sigma\} = \sigma \end{aligned}\]

  • Note: The jump of a scalar is vector valued (in the normal direction), the jump of a vector is scalar

The Primal Formulation

  • Summing over all \(K\), the flux formulation can be written \[\begin{aligned} &\int_\Omega \sigma_h\cdot\tau\,dx = -\int_\Omega u_h \nabla_h\cdot\tau\,dx + \int_\Gamma \left[\left[\right.\right.\hat{u} \left.\left.\right]\right]\cdot \{\tau\}\,ds +\int_{\Gamma^0} \{\hat{u}\} \left[\left[\right.\right.\tau\left.\left.\right]\right]\,ds \\ &\int_\Omega \sigma_h\cdot\nabla_h v\,dx - \int_\Gamma \{\hat{\sigma}\}\cdot\left[\left[\right.\right.v \left.\left.\right]\right]\,ds - \int_{\Gamma^0} \left[\left[\right.\right.\hat{\sigma} \left.\left.\right]\right]\{v\} = \int_\Omega f\,v\,dx \end{aligned}\]

  • With some manipulations, \(\sigma_h\) can be expressed as \[\begin{aligned} \sigma_h = \sigma_h(u_h) := \nabla_h u_h -r(\left[\left[\right.\right.\hat{u}(u_h)-u_h \left.\left.\right]\right]) - l(\{\hat{u}(u_h)-u_h\}) \end{aligned}\] where \(r,l\) are lifting operators defined by \[\begin{aligned} \hspace{-5mm} \int_\Omega r(\phi)\cdot\tau\,dx = -\int_\Gamma \phi\cdot\{\tau\}\,ds, \quad \int_\Omega l(q)\cdot\tau\,dx = -\int_{\Gamma^0} q\left[\left[\right.\right.\tau \left.\left.\right]\right]\,ds \quad \forall \tau\in\Sigma_h \end{aligned}\]

The Primal Formulation

  • This leads to the primal formulation \[\begin{aligned} B_h(u_h,v) &= \int_\Omega f\,v\,dx \qquad \forall v\in V_h \end{aligned}\] with the primal form \[\begin{aligned} B_h(u_h,v) &= \int_\Omega \nabla_h u_h \cdot \nabla_h v\,dx +\int_\Gamma (\left[\left[\right.\right.\hat{u}-u_h \left.\left.\right]\right]\cdot \{\nabla_h v\} -\{\hat{\sigma}\}\cdot\left[\left[\right.\right.v\left.\left.\right]\right])\,ds \\ & \qquad \qquad + \int_{\Gamma^0} (\{\hat{u}-u_h\} \left[\left[\right.\right.\nabla_h v\left.\left.\right]\right] -\left[\left[\right.\right.\hat{\sigma}\left.\left.\right]\right]\{v\} )\,ds \end{aligned}\]

  • Standard FEM formulation without \(\sigma_h\)

  • In implementations it is often easier to work directly with the flux formulation

Consistency and Conservation

  • The numerical fluxes are consistent if for smooth functions \(v\) \[\begin{aligned} \hat{u}(v) = v|_\Gamma ,\qquad \hat{\sigma}(v,\nabla v) = \nabla v|_\Gamma \end{aligned}\] \(\Rightarrow\) consistency of the primal formulation and Galerkin orthogonality \(B_h(u-u_h,v)=0\), \(\forall v\in V_h\)

  • The numerical fluxes are conservative if \(\hat{u}(\cdot)\) and \(\hat{\sigma}(\cdot,\cdot)\) are single-valued on \(\Gamma\)
    \(\Rightarrow\) adjoint consistency of the primal form

Some DG Methods

  • Some of the most important schemes are summarized below:

Method \(\hat{u}_K\) \(\hat{\sigma}_K\) Stable
Bassi-Rebay (BR1) \(\{u_h\}\) \(\{\sigma_h\}\) \(\times\)
Bassi-Rebay (BR2) \(\{u_h\}\) \(\{\nabla_h u_h\} - \alpha_r(\left[\left[\right.\right.u_h \left.\left.\right]\right])\) \(\inf_e \eta_e>3\)
Interior Penalty \(\{u_h\}\) \(\{\nabla_h u_h\} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]\) \(C_{11}>C_{11}^*\)
LDG \(\{u_h\} +\boldsymbol{C}_{12} \cdot \left[\left[\right.\right.u_h \left.\left.\right]\right]\) \(\{\sigma_h\} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]- \boldsymbol{C}_{12} \left[\left[\right.\right.\sigma_h \left.\left.\right]\right]\) \(C_{11}>0\)
  • \(\alpha_r(\phi) = -\eta_e \{r_e(\phi)\}\) on an edge \(e\), where \(r_e\) is defined by \[\begin{aligned} \int_\Omega r_e(\varphi) \cdot\tau\,dx = -\int_e \varphi\cdot \{\tau\}\,ds, \qquad \forall\tau \in\Sigma_h, \varphi\in [L^1(e)]^2 \end{aligned}\]

  • \(C_{11}^*\) is mesh dependent, explicit form derived by Shahbazi (2005)

  • The methods BR2, IP, and LDG are all commonly used

The Local Discontinuous Galerkin (LDG) Method

The LDG Method

  • In the LDG method, we use the fluxes \[\begin{aligned} \hat{\sigma}_K &= \{\sigma_h\} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]- \boldsymbol{C}_{12} \left[\left[\right.\right.\sigma_h \left.\left.\right]\right]\\ \hat{u}_K &= \{u_h\} +\boldsymbol{C}_{12} \cdot \left[\left[\right.\right.u_h \left.\left.\right]\right] \end{aligned}\]

  • Here, \(C_{11}>0\) (or zero for the minimal dissipation LDG method, Cockburn and Dong 2007)


  • An important special case for \(\boldsymbol{C}_{12}\) is the choice \[\begin{aligned} \boldsymbol{C}_{12} = \frac12 (S_{K^+}^{K^-1} n^+ + S_{K^-}^{K^+} n^- ) \end{aligned}\] where \(S_{K^+}^{K^-}\in\{0,1\}\) is a switch for the edge shared by \(K^-\) and \(K^+\)

  • This leads to a simple upwind/downwind scheme: \[\begin{aligned} \hat{\sigma}_K = C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]+ \begin{cases} \sigma_h^+ & \text{if }S_{K^+}^{K^-}=0 \\ \sigma_h^- & \text{if }S_{K^+}^{K^-}=1 \end{cases}, \qquad \hat{u}_K = \begin{cases} u_h^- & \text{if }S_{K^+}^{K^-}=0 \\ u_h^+ & \text{if }S_{K^+}^{K^-}=1 \end{cases} \end{aligned}\]

LDG Switch Functions

  • Natural switch: Order the elements, let \(N_K\) be the index of element \(K\), and set \[\begin{aligned} S_{K^+}^{K^-} = 1 \text{ if } N_{K^+}>N_{K^-},\quad 0 \text{ otherwise}. \end{aligned}\] Simple, leads to beneficial matrix structure, but unstable if \(C_{11}=0\) in the original LDG method

  • Consistent switch: For example define any constant vector \(\beta\) and set \[\begin{aligned} S_{K^+}^{K^-} = 1 \text{ if } n^+ \cdot \beta>0,\quad 0 \text{ otherwise}. \end{aligned}\]

  • In general, any choice of switch leads to a non-compact stencil


Natural switch


Consistent switch

The Compact Discontinuous Galerkin (CDG) Method

The Compact DG (CDG) Method

  • To address the non-compactness of the LDG method and its sensitivity to the switch, Peraire and Persson developed the Compact DG method (2008)

  • Recall the original LDG fluxes: \[\begin{aligned} \hat{\sigma}_K &= \{\sigma_h\} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]- \boldsymbol{C}_{12} \left[\left[\right.\right.\sigma_h \left.\left.\right]\right]\\ \hat{u}_K &= \{u_h\} +\boldsymbol{C}_{12} \cdot \left[\left[\right.\right.u_h \left.\left.\right]\right] \end{aligned}\]

  • Now, introduce the edge fluxes \(\sigma_h^e\) on edge \(e\) by \[\begin{aligned} \int_K \sigma_h^e\cdot\tau\,dx &= -\int_K u_h \nabla\cdot\tau\,dx + \int_{\partial K} \hat{u}^e_K\,n_K\cdot\tau\,ds, && \forall\tau\in[\mathcal{P}_p(K)]^2 \end{aligned}\]

  • where \[\begin{aligned} \hat{u}^e_K = \begin{cases} \hat{u}_K & \text{on edge }e\text{, as defined above} \\ u_h & \text{otherwise} \end{cases} \end{aligned}\]

The Compact DG (CDG) Method

  • The numerical fluxes for CDG are then simply given by \[\begin{aligned} \hat{\sigma}^e_K &= \{\sigma_h^e\} + C_{11} \left[\left[\right.\right.u_h \left.\left.\right]\right]- \boldsymbol{C}_{12} \left[\left[\right.\right.\sigma_h^e \left.\left.\right]\right] \end{aligned}\] on edge \(e\)

  • The modification eliminates the non-compact terms in the primal form, while retaining all the good properties of the LDG scheme

  • In addition, better stability properties are observed with in particular less sensitivity to the choice of switch function  


The CDG Method – Summary




  • Element-wise compact stencil

  • Less connectivities than LDG/BR2/IP

  • More accurate than LDG and BR2


Properties of the CDG Method

Switches and Null-space Dimensions

  • Unlike the LDG scheme, the CDG scheme appears to be stable for \(C_{11}=0\) and an inconsistent switch such as highest element number

  • Simple test [Sherwin et al 05]: Poisson problem, periodic boundary conditions, expected nullspace dimension = 1


Nullspace dimension

Polynomial order \(p\) \(1\) \(2\) \(3\) \(4\) \(5\) \(6\) \(7\)
Consistent switch CDG \(1\) \(1\) \(1\) \(1\) \(1\) \(1\) \(1\)
LDG \(1\) \(1\) \(1\) \(1\) \(1\) \(1\) \(1\)
Natural switch CDG \(1\) \(1\) \(1\) \(1\) \(1\) \(1\) \(1\)
LDG \(3\) \(4\) \(5\) \(6\) \(7\) \(8\) \(9\)

ILU and Switch Orientation

  • Orientation of lower-triangular blocks important for ILU sparsity

  • Take advantage of CDG’s insensitivity to orientation

image \(A\)


image \(L\)

image \(U\)

Switch 1:
Same LU storage

image \(A\)


image \(L\)

image \(U\)

Switch 2:
More LU storage

Switches and Null-space Dimensions

  • No additional non-zeros in block-ILU(0) factorization using CDG

  • Dense lower-triangular blocks using BR2 / IP



Stiffness Matrix

640 non-zeros

Block ILU(0)

640 non-zeros

Switches and Null-space Dimensions

  • No additional non-zeros in block-ILU(0) factorization using CDG

  • Dense lower-triangular blocks using BR2 / IP

BR2 / IP


Stiffness Matrix

784 non-zeros

Block ILU(0)

892 non-zeros

Matrix Representation

  • Block matrix representation fundamental for high performance

    • Solver algorithms based on blocks

    • Up to 10 times higher performance with optimized BLAS

  • Compact stencil \(\Longrightarrow\) Matrix structure given by mesh connectivities

  • Hard to store LDG/BR2/IP efficiently

CDG – 2 arrays

LDG – 3 arrays + struct

BR2 / IP – 3 arrays