# Discontinuous Galerkin Methods for Conservation Laws

Math 228B Numerical Solutions of Differential Equations

The Finite Volume Method = Galerkin FEM

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

• Look for solutions in space of piecewise constant functions $$V_h$$ \begin{aligned} u_h(x) = \sum_{k=1}^n u_k \varphi_k(x),\qquad \varphi_k(x) = \left\{ \begin{array}{rl} 1 & x_{k-1}<x<x_k \\ 0 & \text{otherwise} \end{array} \right. \end{aligned}

The Finite Volume Method = Galerkin FEM

• Galerkin formulation: Find $$u_h\in V_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,\quad \forall v\in V_h \end{aligned}

• Set $$v=\varphi_k= \begin{cases} 1 & x\in[x_{k-1},x_k] \\ 0 & \text{otherwise} \end{cases}$$ \begin{aligned} \hspace{-1.6cm} \int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t}\,dx +\int_{x_{k-1}}^{x_k} \frac{\partial f(u_h)}{\partial x}\,dx = 0 \Longleftrightarrow \int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t}\,dx +\left[ f(u_h(x))\right]_{x_{k-1}}^{x_k} = 0 \end{aligned}

• Since $$u_h$$ is discontinuous at $$x_k$$ and $$x_{k-1}$$, use a numerical flux function $$F(u_R,u_L)$$ to obtain: \begin{aligned} h\frac{\partial u_k}{\partial t} + F(u_{k+1},u_{k}) - F(u_k,u_{k-1}) = 0 \end{aligned}

• This is a standard finite volume method on a uniform grid

The Discontinuous Galerkin Method

• Generalize the Galerkin FEM approach to the space of piecewise polynomials of degree $$p$$

• 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 \varphi_i^k(x) \end{aligned}

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

The Discontinuous Galerkin Method

• Galerkin formulation: Find $$u_h\in V_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,\quad \forall v\in V_h \end{aligned}

• Set $$v=\varphi_i^k$$ and integrate by parts \begin{aligned} \hspace{-10mm} \int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t} \varphi_i^k\,dx + \left[ f(u_h(x))\varphi_i^k(x) \right]_{x_{k-1}}^{x_k} -\int_{x_{k-1}}^{x_k} f(u_h(x))\frac{d\varphi_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} \varphi_i^k\,dx + F(u_0^{k+1},u_p^k)\varphi_i^k(x_{k}) - F(u_0^k,u_p^{k-1})\varphi_i^k(x_{k-1}) \\ &-\int_{x_{k-1}}^{x_k} f(u_h(x))\frac{d\varphi_i^k}{dx}\,dx = 0 \end{aligned}

The Discontinuous Galerkin Method

• 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_{j=0}^p u_j^k \varphi_j^k(x) \right) \varphi_i^k\,dx - \int_{x_{k-1}}^{x_k} \left( \sum_{j=0}^p u_j^k \varphi_j^k(x) \right) \frac{d\varphi_i^k}{dx}\,dx \\ & +u_p^k\varphi_i^k(x_k)-u_p^{k-1}\varphi_i^k(x_{k-1}) = 0 \end{aligned}

• Rearrange to obtain a linear system of equations \begin{aligned} M^k \dot{\boldsymbol{u}}^k-C^k\boldsymbol{u}^k + \begin{pmatrix} -u_p^{k-1} \\ 0 \\ \vdots \\ 0 \\ u_p^k \end{pmatrix} = 0 \end{aligned} for element $$k$$, with elementary matrices \begin{aligned} M_{ij}^k = \int_{x_{k-1}}^{x_k} \varphi_i^k\varphi_j^k\,dx\text{ and } C_{ij}^k = \int_{x_{k-1}}^{x_k} \frac{d\varphi_i^k}{dx}\varphi_j^k\,dx \end{aligned}

Calculating Elementary Matrices

• Consider an element of degree $$p$$, width $$h$$, and a nodal basis for the points $$s_i$$, $$i=0,\ldots,p$$

• Equidistant points $$s_i=ih/p$$ only good for low $$p$$

• Better choice: Chebyshev or Gauss-Lobatto nodes

• Write basis functions as $$\varphi_i(s) = \sum_{j=0}^p c_i^j P_j(s)$$, where $$P_j$$ is a basis for the polynomials of degree $$p$$

• Monomial basis $$P_j(s) = s^j$$ only good for low $$p$$

• Better choice: Orthogonal polynomials, e.g. Legendre

• Nodal basis functions are defined by \begin{aligned} \varphi_i(s_j) = \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} \hspace{-12mm} \begin{pmatrix} P_0(s_0) & P_1(s_0) & \cdots & P_p(s_0) \\ P_0(s_1) & P_1(s_1) & \cdots & P_p(s_1) \\ \vdots & \vdots & \ddots & \vdots \\ P_0(s_p) & P_1(s_p) & \cdots & P_p(s_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 \varphi_i(s)\varphi_j(s)\, ds \\ C_{ij}&=\int_0^h \varphi_i'(s)\varphi_j(s)\, ds \end{aligned}

The DG method – General systems of conservation laws

• (Reed/Hill 1973, Lesaint/Raviart 1974, Cockburn/Shu 1989-)

• Consider a first-order system of conservation laws: \begin{aligned} \boldsymbol{u}_t + \nabla \cdot \boldsymbol{F}(\boldsymbol{u}) = 0 \end{aligned}

• Triangulate domain $$\Omega$$ into elements $$\kappa\in \mathcal{T}_h$$

• Seek approximate solution $$\boldsymbol{u}_h$$ in space of element-wise polynomials: \begin{aligned} \boldsymbol{V}_h^p = \{\boldsymbol{v}\in L^2(\Omega) : \boldsymbol{v}|_\kappa\in P^p(\kappa) \,\,\forall \kappa \in T_h \} \end{aligned}

• Multiply by test function $$\boldsymbol{v}_h\in\boldsymbol{V}_h^p$$, integrate over element $$\kappa$$: \begin{aligned} \int_\kappa \left[ (\boldsymbol{u}_h)_t + \nabla \cdot \boldsymbol{F}(\boldsymbol{u}_h) \right] \boldsymbol{v}_h\, d\boldsymbol{x} = 0 \end{aligned}

The DG method – General systems of conservation laws

• Integrate by parts: \begin{aligned} \hspace{-5mm} \int_\kappa \left[(\boldsymbol{u}_h)_t\right] \boldsymbol{v}_h\,d\boldsymbol{x} - \int_\kappa \boldsymbol{F}(\boldsymbol{u}_h) \nabla\boldsymbol{v}_h\,d\boldsymbol{x} + \int_{\partial\kappa} \hat{\boldsymbol{F}}(\boldsymbol{u}_h^+,\boldsymbol{u}_h^-,\hat{\boldsymbol{n}}) \boldsymbol{v}_h^+\,ds = 0 \end{aligned} with numerical flux function $$\hat{\boldsymbol{F}}(\boldsymbol{u}_L,\boldsymbol{u}_R,\hat{\boldsymbol{n}})$$ for left/right states $$\boldsymbol{u}_L$$,$$\boldsymbol{u}_R$$ in direction $$\hat{\boldsymbol{n}}$$ (Godunov, Roe, Osher, Van Leer, Lax-Friedrichs, etc)

• Global problem: Find $$\boldsymbol{u}_h\in\boldsymbol{V}_h^p$$ such that this weighted residual is zero for all $$\boldsymbol{v}_h\in\boldsymbol{V}_h^p$$

• Error $$= \mathcal{O}(h^{p+1})$$ for smooth solutions

The DG Method – Observations

• Reduces to the finite volume method for $$p=0$$: \begin{aligned} (\boldsymbol{u}_h)_t A_\kappa + \int_{\partial\kappa} \hat{\boldsymbol{F}}(\boldsymbol{u}_h^+, \boldsymbol{u}_h^-,\hat{\boldsymbol{n}})\,ds = 0 \end{aligned}

• Boundary conditions enforced naturally for any degree $$p$$

• Block-diagonal mass matrix (no overlap between basis functions)

• Block-wise compact stencil – neighboring elements connected

Mass Matrix

Jacobian

Convection-Diffusion, the LDG method

• Consider the convection-diffusion equation \begin{aligned} \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} - \mu\frac{\partial ^2u}{\partial x^2} = 0 \end{aligned}

• Split into system of first order equations: \begin{aligned} \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} - \mu\frac{\partial \sigma}{\partial x} &= 0 \\ \frac{\partial u}{\partial x} &= \sigma \end{aligned}

• Galerkin formulation: Find $$u_h,\sigma_h \in V_h$$ such that \begin{aligned} &\int_0^1 \frac{\partial u_h}{\partial t}v\, dx + \int_0^1 \left(\frac{\partial f(u_h)}{\partial x} - \mu\frac{\partial \sigma_h}{\partial x} \right) v\,dx = 0,\quad \forall v\in V_h \\ &\int_0^1 \frac{\partial u_h}{\partial x}\tau\,dx = \int_0^1 \sigma_h \tau\,dx,\quad \forall \tau\in V_h \end{aligned}

Convection-Diffusion, the LDG method

• Set $$v,\tau=\varphi_i^k$$ and integrate by parts \begin{aligned} \hspace{-10mm} &\int_{x_{k-1}}^{x_k} \frac{\partial u_h}{\partial t} \varphi_i^k\,dx + \left[ \left(f(u_h(x)) - \mu\sigma_h(x)\right)\varphi_i^k(x) \right]_{x_{k-1}}^{x_k} \\ & \qquad\qquad-\int_{x_{k-1}}^{x_k} \left(f(u_h(x))-\mu\sigma_h(x) \right)\frac{d\varphi_i^k}{dx}\,dx = 0, \qquad \forall i,k\\ &\left[ u_h(x)\varphi_i^k(x) \right]_{x_{k-1}}^{x_k} \\ & \qquad\qquad -\int_{x_{k-1}}^{x_k} u_h(x) \frac{d\varphi_i^k}{dx}\,dx = \int_{x_{k-1}}^{x_k} \sigma_h(x) \varphi_i^k \,dx,\qquad \forall i,k \end{aligned}

• Use numerical flux functions $$\hat{f}(u_R,u_L)$$, $$\hat{\sigma}(\sigma_R,\sigma_L)$$, $$\hat{u}(u_R,u_L)$$ at the discontinuities

• Example: $$f(u)=u$$, $$\hat{f}(u_R,u_L) = u_L$$, $$\hat{\sigma}(\sigma_R,\sigma_L) = \sigma_L$$, $$\hat{u}(u_R,u_L) = u_R$$ (upwinding for the convection, LDG upwinding/downwinding for the diffusion)

Convection-Diffusion, the LDG method

• After discretization, this leads to the ODEs \begin{aligned} &M^k \dot{\boldsymbol{u}}^k-C^k\left(\boldsymbol{u}^k - \mu\boldsymbol{\sigma}^k \right) + \begin{pmatrix} -u_p^{k-1} + \mu \sigma_p^{k-1} \\ 0 \\ \vdots \\ 0 \\ u_p^k - \mu \sigma_p^k \end{pmatrix} = 0 \\ & M^k \boldsymbol{\sigma}^k = -C^k \boldsymbol{u}^k + \begin{pmatrix} -u_0^{k} \\ 0 \\ \vdots \\ 0 \\ u_0^{k+1} \end{pmatrix} \end{aligned}

• For each element $$k$$, first solve for $$\boldsymbol{\sigma}^k$$, then insert into main equation as before