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\)):
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
\(A\)
\(=\)
\(L\)
\(U\)
Switch 1:
Same LU storage
\(A\)
\(=\)
\(L\)
\(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
CDG
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