Discontinuous Galerkin method

Much like the continuous Galerkin (CG) method, the discontinuous Galerkin (DG) method is a finite element method formulated relative to a weak formulation of a particular model system. Unlike traditional CG methods that are conforming, the DG method works over a trial space of functions that are only piecewise continuous, and thus often comprise more inclusive function spaces than the finite-dimensional inner product subspaces utilized in conforming methods. As an example, consider the continuity equation for a scalar unknown \rho in a spatial domain \Omega without "sources" or "sinks" : : \frac{\partial \rho} {\partial t} + \nabla \cdot \mathbf{J} = 0, where \mathbf{J} is the flux of \rho. Now consider the finite-dimensional space of discontinuous piecewise polynomial functions over the spatial domain \Omega restricted to a discrete triangulation \Omega_h , written as : S_h^p(\Omega_h)=\{v_{|\Omega_{e_i}}\in P^p(\Omega_{e_i}), \ \ \forall\Omega_{e_i}\in \Omega_h\} for P^p(\Omega_{e_i}) the space of polynomials with degrees less than or equal to p over element \Omega_{e_i} indexed by i. Then for finite element shape functions N_j\in P^p the solution is represented by :\rho_h^i = \sum_{j=1}^{\text{dofs}} \rho_j^i (t)N_j^i (\boldsymbol{x}), \quad \forall \boldsymbol{x} \in \Omega_{e_i}. Then similarly choosing a test function :\varphi_h^i(\boldsymbol{x})=\sum_{j=1}^{\text{dofs}} \varphi_j^i N_j^i(\boldsymbol{x}), \quad \forall \boldsymbol{x}\in\Omega_{e_i}, multiplying the continuity equation by \varphi_h^i and integrating by parts in space, the semidiscrete DG formulation becomes: : \frac{d}{dt}\int_{\Omega_{e_i}}\rho_h^i\varphi_h^i \, d\boldsymbol{x} + \int_{\partial\Omega_{e_{i}}} \varphi_h^i \mathbf{J}_h \cdot\boldsymbol{n} \, d\boldsymbol{x} = \int_{\Omega_{e_i}}\mathbf{J}_h\cdot\nabla\varphi_h^i \, d\boldsymbol{x}. == Scalar hyperbolic conservation law ==

A scalar hyperbolic conservation law is of the form : \begin{align} \partial_t u + \partial_x f(u) &= 0 \quad \text{for} \quad t>0,\, x\in \R \\ u(0,x) &= u_0(x)\,, \end{align} where one tries to solve for the unknown scalar function u \equiv u(t,x) , and the functions f,u_0 are typically given. Space discretization The x -space will be discretized as : \R = \bigcup_k I_k \,, \quad I_k := \left( x_k, x_{k+1} \right) \quad \text{for} \quad x_k Furthermore, we need the following definitions : h_k := | I_k | \,, \quad h := \sup_k h_k \,, \quad \hat{x}_k := x_k + \frac{h_k}{2}\,. Basis for function space We derive the basis representation for the function space of our solution u . The function space is defined as : S_h^p := \left\lbrace v \in L^2(\R) : v\Big|_{I_k} \in \Pi_p \right\rbrace \quad \text{for} \quad p \in \N_0 \,, where {v|}_{I_k} denotes the restriction of v onto the interval I_k , and \Pi_p denotes the space of polynomials of maximal degree p . The index h should show the relation to an underlying discretization given by \left(x_k\right)_k . Note here that v is not uniquely defined at the intersection points (x_k)_k . At first we make use of a specific polynomial basis on the interval [-1,1] , the Legendre polynomials (P_n)_{n\in\N_0} , i.e., : P_0(x) = 1 \,,\quad P_1(x)=x \,,\quad P_2(x) = \frac{1}{2} (3x^2-1) \,,\quad \dots Note especially the orthogonality relations : \left\langle P_i,P_j \right\rangle_{L^2([-1,1])} = \frac{2}{2i+1} \delta_{ij} \quad \forall \, i,j \in \N_0 \,. Transformation onto the interval [0,1] , and normalization is achieved by functions (\varphi_i)_i : \varphi_i (x) := \sqrt{2i+1} P_i(2x-1) \quad \text{for} \quad x\in [0,1]\,, which fulfill the orthonormality relation : \left\langle \varphi_i,\varphi_j \right\rangle_{L^2([0,1])} = \delta_{ij} \quad \forall \, i,j \in \N_0 \,. Transformation onto an interval I_k is given by \left( \bar{\varphi}_{ki}\right)_i : \bar{\varphi}_{ki} := \frac{1}{\sqrt{h_k}} \varphi_i \left( \frac{x-x_k}{h_k} \right) \quad \text{for} \quad x\in I_k\,, which fulfill : \left\langle \bar{\varphi}_{ki},\bar{\varphi}_{kj} \right\rangle_{L^2(I_k)} = \delta_{ij} \quad \forall \, i,j \in \N_0 \forall \, k \,. For L^\infty -normalization we define \varphi_{ki}:= \sqrt{h_k} \bar{\varphi}_{ki} , and for L^1 -normalization we define \tilde{\varphi}_{ki}:= \frac{1}{\sqrt{h_k}} \bar{\varphi}_{ki} , s.t. : \| \varphi_{ki} \|_{L^\infty (I_k) } = \| \varphi_i \|_{L^\infty ([0,1]) } =: c_{i,\infty} \quad \text{and} \quad \| \tilde{\varphi}_{ki} \|_{L^1 (I_k) } = \| \varphi_i \|_{L^1 ([0,1]) } =: c_{i,1} \,. Finally, we can define the basis representation of our solutions u_h : \begin{align} u_h(t,x) :=& \sum_{i=0}^p u_{ki}(t) \varphi_{ki} (x) \quad \text{for} \quad x \in (x_k,x_{k+1}) \\ u_{ki} (t) =& \left\langle u_h(t, \cdot ),\tilde{\varphi}_{ki} \right\rangle_{L^2(I_k)} \,. \end{align} Note here, that u_h is not defined at the interface positions. Besides, prism bases are employed for planar-like structures, and are capable for 2-D/3-D hybridation. DG-scheme The conservation law is transformed into its weak form by multiplying with test functions, and integration over test intervals : \begin{align} \partial_t u + \partial_x f(u) &= 0 \\ \Rightarrow \quad \left\langle \partial_t u , v \right\rangle_{L^2(I_k)} + \left\langle \partial_x f(u) , v \right\rangle_{L^2(I_k)} &= 0 \quad \text{for} \quad \forall \, v \in S_h^p \\ \Leftrightarrow \quad \left\langle \partial_t u , \tilde{\varphi}_{ki} \right\rangle_{L^2(I_k)} + \left\langle \partial_x f(u) , \tilde{\varphi}_{ki} \right\rangle_{L^2(I_k)} &= 0 \quad \text{for} \quad \forall \, k \; \forall\, i \leq p \,. \end{align} By using partial integration one is left with : \begin{align} \frac{\mathrm d}{\mathrm d t} u_{ki}(t) + f(u(t, x_{k+1} )) \tilde{\varphi}_{ki}(x_{k+1}) - f(u(t, x_k )) \tilde{\varphi}_{ki}(x_k) - \left\langle f(u(t,\,\cdot\,)) , \tilde{\varphi}_{ki}' \right\rangle_{L^2(I_k)} =0 \quad \text{for} \quad \forall \, k \; \forall\, i \leq p \,. \end{align} The fluxes at the interfaces are approximated by numerical fluxes g with : g_k := g(u_k^-,u_k^+) \,, \quad u_k^\pm := u(t,x_k^\pm) \,, where u_k^{\pm} denotes the left- and right-hand sided limits. Finally, the DG-Scheme can be written as : \begin{align} \frac{\mathrm d}{\mathrm d t} u_{ki}(t) + g_{k+1} \tilde{\varphi}_{ki}(x_{k+1}) - g_k \tilde{\varphi}_{ki}(x_k) - \left\langle f(u(t,\,\cdot\,)) , \tilde{\varphi}_{ki}' \right\rangle_{L^2(I_k)} =0 \quad \text{for} \quad \forall \, k \; \forall\, i \leq p \,. \end{align} == Scalar elliptic equation ==

A scalar elliptic equation is of the form : \begin{align} -\partial_{xx} u &= f(x) \quad \text{for} \quad x\in (a,b) \\ u(x) &= g(x)\,\quad\text{for}\,\quad x=a,b \end{align} This equation is the steady-state heat equation, where u is the temperature. Space discretization is the same as above. We recall that the interval (a,b) is partitioned into N+1 intervals of length h. We introduce jump [{}\cdot{}] and average \{{}\cdot{}\} of functions at the node x_k: : [v]\Big|_{x_k} = v(x_k^+)-v(x_k^-), \quad \{v\}\Big|_{x_k} = 0.5 (v(x_k^+)+v(x_k^-)) The interior penalty discontinuous Galerkin (IPDG) method is: find u_h satisfying : A(u_h,v_h) + A_{\partial}(u_h,v_h) = \ell(v_h) + \ell_\partial(v_h) where the bilinear forms A and A_\partial are : A(u_h,v_h) = \sum_{k=1}^{N+1} \int_{x_{k-1}}^{x_k}\partial_x u_h \partial_x v_h -\sum_{k=1}^N \{ \partial_x u_h\}_{x_k} [v_h]_{x_k} +\varepsilon\sum_{k=1}^N \{ \partial_x v_h\}_{x_k} [u_h]_{x_k} +\frac{\sigma}{h} \sum_{k=1}^N [u_h]_{x_k} [v_h]_{x_k} and : A_\partial(u_h,v_h) = \partial_x u_h(a) v_h(a) -\partial_x u_h(b) v_h(b) -\varepsilon \partial_x v_h(a) u_h(a) + \varepsilon\partial_x v_h(b) u_h(b) +\frac{\sigma}{h} \big(u_h(a) v_h(a) + u_h(b) v_h(b)\big) The linear forms \ell and \ell_\partial are : \ell(v_h) = \int_a^b f v_h and : \ell_\partial(v_h) = -\varepsilon \partial_x v_h(a) g(a) + \varepsilon\partial_x v_h(b) g(b) +\frac{\sigma}{h} \big( g(a) v_h(a) + g(b) v_h(b) \big) The penalty parameter \sigma is a positive constant. Increasing its value will reduce the jumps in the discontinuous solution. The term \varepsilon is chosen to be equal to -1 for the symmetric interior penalty Galerkin method; it is equal to +1 for the non-symmetric interior penalty Galerkin method. ==Direct discontinuous Galerkin method==

The direct discontinuous Galerkin (DDG) method is a new discontinuous Galerkin method for solving diffusion problems. In 2009, Liu and Yan first proposed the DDG method for solving diffusion equations. The advantage of this method compared with Discontinuous Galerkin method is that the direct discontinuous Galerkin method derives the numerical format by directly taking the numerical flux of the function and the first derivative term without introducing intermediate variables. We can still get reasonable numerical results by using this method, and due to the simpler derivation process, the amount of calculation is greatly reduced. The direct discontinuous finite element method is a branch of the Discontinuous Galerkin methods. It mainly includes transforming the problem into variational form, regional unit splitting, constructing basis functions, forming and solving discontinuous finite element equations, and convergence and error analysis. For example, consider a nonlinear diffusion equation, which is one-dimensional: : U_t - {(a(U)\cdot U_x)}_x = 0 \ \ in \ (0,1) \times (0,T), in which U(x,0) = U_0(x)\ \ on\ (0,1) Space discretization Firstly, define \left \{ I_j =\left ( x_{j-\frac{1}{2}},\ x_{j+\frac{1}{2}} \right ), j=1...N\right \}, and \Delta x_j=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}. Therefore we have done the space discretization of x. Also, define \Delta x=\max_{1\leq j. We want to find an approximation u to U such that \forall t\in \left [ 0,T \right ], u \in \mathbb{V}_{\Delta x}, \mathbb{V}_{\Delta x} := \left \{ v\in L^2\left ( 0,1 \right ):{v|}_{I_j}\in P^k\left ( I_j \right ), \ j=1,...,N \right \}, P^k\left ( I_j \right ) is the polynomials space in I_j with degree at most k. Formulation of the scheme Flux: h:= h\left ( U,U_x \right )=a\left ( U \right )U_x. U : the exact solution of the equation. Multiply the equation with a smooth function v\in H^1\left (0,1 \right ) so that we obtain the following equations: \int _{I_j} U_tvdx-h_{j+\frac{1}{2}}v_{j+\frac{1}{2}}+h_{j-\frac{1}{2}} v_{j-\frac{1}{2}} +\int a\left ( U \right )U_xv_xdx=0 , \int _{I_j} U\left ( x,0 \right )v\left ( x \right )dx=\int _{I_j}U_0\left ( x \right )v\left ( x \right )dx Here v is arbitrary, the exact solution U of the equation is replaced by the approximate solution u, that is to say, the numerical solution we need is obtained by solving the differential equations. The numerical flux Choosing a proper numerical flux is critical for the accuracy of DDG method. The numerical flux needs to satisfy the following conditions: ♦ It is consistent with h={b\left ( u \right )}_x=a\left ( u \right )u_x ♦ The numerical flux is conservative in the single value on x_{j+\frac{1}{2}}. ♦ It has the L^2-stability; ♦ It can improve the accuracy of the method. Thus, a general scheme for numerical flux is given: :\widehat{h}=D_xb(u)=\beta_0\frac{\left [ b\left ( u \right ) \right ]}{\Delta x}+\overline{{b\left ( u \right )}_x}+\sum_{m=1}^{\frac{k}{2}}\beta_m{\left ( \Delta x \right )}^{2m-1}\left [ \partial _x^{2m}b\left ( u \right ) \right ] In this flux, k is the maximum order of polynomials in two neighboring computing units. \left [\cdot \right ] is the jump of a function. Note that in non-uniform grids, \Delta x should be \left ( \frac{\Delta x_j+\Delta x_{j+1}}{2} \right ) and \frac{1}{N} in uniform grids. Error estimates Denote that the error between the exact solution U and the numerical solution u is e= u-U . We measure the error with the following norm: \left | \left | \left | v(\cdot ,t) \right | \right | \right | ={\left (\int_{0}^{1}v^2dx+\left ( 1-\gamma \right )\int_{0}^{t}\sum_{j=1}^{N}\int_{I_j} v_x^2 dxd\tau +\alpha \int_{0}^{t}\sum_{j=1}^{N} {\left [ v \right ] }^2/\Delta x\cdot d\tau \right )}^{0.5} and we have \left | \left | \left | U(\cdot ,T) \right | \right | \right |\leq \left | \left | \left | U(\cdot ,0) \right | \right | \right |,\left | \left | \left | u(\cdot ,T) \right | \right | \right |\leq \left | \left | \left | U(\cdot ,0) \right | \right | \right | ==See also==