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\)):

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 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

• 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:

• \(\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 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

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

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

ILU and Switch Orientation

• Orientation of lower-triangular blocks important for ILU sparsity

• Take advantage of CDG’s insensitivity to orientation

Switches and Null-space Dimensions

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

• Dense lower-triangular blocks using BR2 / IP

Switches and Null-space Dimensions

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

• Dense lower-triangular blocks using BR2 / IP

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