% beamer version of chapter 7 "The Discontinuous GalerkinMethod" in "Computational Seismology, A Practical Introduction" (H. Igel)

\documentclass[aspectratio=169]{beamer}

\mode<presentation>
{
  \usetheme{metropolis}
	\usecolortheme{crane}

  \setbeamercovered{transparent}
  \setbeamertemplate{navigation symbols}{}
}


\usepackage[english]{babel}
\usepackage[latin1]{inputenc}

\usepackage{times}
\usepackage[T1]{fontenc}
% Or whatever. Note that the encoding and the font should match. If T1
% does not look nice, try deleting the line with the fontenc.

% Math packages
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{graphicx}
%using a different letter style for capital letters
\usepackage{mathrsfs}

\usepackage{mdframed}
\usepackage{listings}
\usepackage[percent]{overpic}
\usepackage{capt-of}
\usepackage{tikz}
\usepackage{grffile}

% Own definitions ... 
\newcommand{\eqdef}{\mathrel{\mathop:}=} % for symbol :=
\newcommand*\Eval[3]{\left.#1\right\rvert_{#2}^{#3}} % for integration bound |
\newcommand{\figpath}{c:/mydata/Heiner/Text/Books/OUP/TheBook/Part2/Figs_dg}
\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Discontinuous Galerkin method}

%\subtitle{Introduction} 

\author[Heiner Igel]{Heiner Igel}


\institute[] 
  {Department of Earth and Environmental Sciences \\ Ludwig-Maximilians-University Munich} 



\date[] % (optional)
{}

\subject{Talks}


%\beamerdefaultoverlayspecification{<+->}

\begin{document}

% ----------------------------------------------------------

\begin{frame}
  \titlepage
\end{frame}

% ----------------------------------------------------------

\begin{frame}{Outline}
\fs
  \tableofcontents
  % You might wish to add the option [pausesections]
\end{frame}

% ---------------------------------------------------------

\section{Introduction}

% ---------------------------------------------------------
\subsection{Motivation}

	\begin{frame}{Motivation}
	\begin{minipage}{0.4\textwidth}
		\begin{figure}
			\includegraphics[width=0.8\textwidth]{\figpath/fig_dg_grid.png}
		\end{figure}
	\end{minipage}
	\mbox{}		
	\hfill		
	\begin{minipage}{0.55\textwidth}
		\begin{itemize}
			\item Efficiently solving the elastic wave equation on tetrahedral grids
			\item Easy implementation of local time stepping
			\item High accuracy of frictional boundary conditions for dynamic rupture problems
		\end{itemize}
	\end{minipage}
	\end{frame}

% ---------------------------------------------------------

\subsection{History}

	\begin{frame}{History}
	\scriptsize
		\begin{itemize}[<+->]
			\item Developed in the Los Alamos National Laboratories for the problem of neutron transport by Reed 1973 formulated on triangular meshes
			\item In the late eighties Cockburn connected discontinuous Galerkin method with high-order Runge-Kutta-type time integration schemes (Cockburn et al., 2000)
			\item First application to elastic wave propagation was published by Kaeser and Dumbser, 2006
			\item Kaeser et al., 2008 carried out a detailed analysis of the convergence properties of the discontinuous Galerkin method
			\item The local time-stepping approach by Dumbser et al., 2007b circumvents the problem of oversampling large parts of the model and thereby reducing the overall computations.
			\item De la Puente et al., 2009b reported a first analysis of scaling and synchronization of the discontinuous Galerkin method which raised the interest of the computational science community in Munich to optimize it.
		\end{itemize}
	\end{frame}
	
% ---------------------------------------------------------

\subsection{The Discontinuous Galerkin Method in a Nutshell}

	\begin{frame}{The Discontinuous Galerkin Method in a Nutshell}
		All concepts that have been discussed previously enter in this method
		\begin{itemize}
			\item[1] Finite-difference extrapolation using e.g. the Euler scheme
			\item[2] Calculation of element-based stiffness and mass matrices
			\item[3] Flux calculations at the element boundaries as encountered in the finite-volume method
			\item[4] Exact nodal interpolation based on Lagrange polynomials
			\item[5] Numerical integration schemes using collocation points
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{Wave Equation}
		\begin{block}{1st order wave equation}
			\begin{equation*}
				\begin{split}
					\rho \partial_t v \ &= \ \partial_x \sigma + f \\
					\partial_t {\sigma} \ &= \ \mu \partial_x v 
				\end{split}
			\end{equation*}
			\end{block}
			\begin{block}{Matrix-vector form}
			\begin{equation*}
				\partial_t Q + A \partial_x Q = f 
			\end{equation*}			
			\end{block}
		\vspace{0.5cm}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{Wave Equation}
		where $\sigma = \sigma_{xy} =  \sigma_{yx}$ representing the only non-zero stress component, and implicitly assuming space-time dependencies, $Q=(v, \sigma)$ is the vector of unknowns and $A$ contains the coefficients of the equation given by
		\begin{equation*}
			A\ = \
				\begin{pmatrix}
					0 & -1/\rho \\
					-\mu & 0 
				\end{pmatrix} \ . 
		\end{equation*}
		$\Longrightarrow$ Linear hyperbolic system, having the same form as the classic advection equation
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method in a Nutshell}
		\begin{itemize}
			\item The Ansatz is to multiply the equation by an arbitrary test function combined with describing the unknown fields with the same set of basis functions (\textbf{Galerkin} principle). 
			\item At the element boundaries the unknown fields can be \textbf{discontinuous}.
			\item Communication between the elements is achieved through a flux scheme based on solutions of the Riemann problem.
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method in a Nutshell}
		\begin{minipage}{0.3\textwidth}
			\begin{figure}
				\includegraphics[width=0.75\textwidth]{\figpath/fig_dg_nutshell.png}
			\end{figure}
		\end{minipage}
	\mbox{}		
	\hfill			
		\begin{minipage}{0.55\textwidth}
		\fs
The wavefield inside each element is described by Lagrange polynomials exactly interpolating at appropriate collocation points. At each time step, a flux term F has to be evaluated at all element boundaries. The extrapolation scheme of the discontinuous Galerkin method can be expressed as 
			\begin{equation*}
				\partial_t Q^{k}(t) \ = \ (M^k)^{-1} (K^k Q^{k}(t) - (F^k_l - F^k_r) Q^{k} (t) )
			\end{equation*}
where $Q(t)$ is the vector of unknowns, $M$ and $K$ are the elemental mass and stiffness matrices, respectively, and $F_{r,l}$ are the flux terms at the left and right boundaries. The upper index $k$ denotes the element (source term is omitted).  
		\end{minipage}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method in a Nutshell}
		Implications of a boundary flux and no global system of equations to solve:
		\begin{itemize}
			\item[1] Obtaining a fully explicit scheme which lends itself to an \textbf{element-based parallelization} scheme
			\item[2] Choice of element size is arbitrary, it has no impact on the solution algorithm (\textbf{h-adaptivity})
			\item[3] The polynomial order in each element can be arbitrarily chosen and no impact on the algorithm (\textbf{p-adaptivity})
			\item[4] Considering the boundary points twice to calculate the fluxes implies an \textbf{increase in the number of degrees of freedom }that of course gets worse with increasing dimensionality
		\end{itemize}
	\end{frame}

% ---------------------------------------------------------

\section{Ingredients}

% ---------------------------------------------------------

	\begin{frame}
		\begin{block}{Source-free version of the coupled elastic wave equation in 1D}
			\begin{equation*}
				\begin{split}
					\partial_t {\sigma} \ - \ \mu \partial_x v \ = \ 0  \\
					\partial_t v - \frac{1}{\rho} \partial_x \sigma \ = \ 0 
				\end{split}
			\end{equation*}
		\end{block}
		\begin{block}{Matrix notation}
			\begin{equation*}
				\partial_t Q + A \partial_x Q = 0 
			\end{equation*}
			where $Q=(\sigma, v)$ is the vector of unknowns and matrix $A$ contains the parameters
			\begin{equation*}
				\textbf{A} = \left(
					\begin{array}{cc}
						0 & -\mu \\
						-\frac{1}{\rho} & 0 \\
					\end{array}
				\right)
			\end{equation*}
		\end{block}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Wave equation as Hyperbolic System}
		In the case of a quadratic matrix $A$ with shape $m \times m$, this leads to an eigenvalue problem. If  we were able to obtain eigenvalues $\lambda_p$ such that
		\begin{equation*}
			\textbf{A} \textbf{x}_p = \lambda_p \textbf{x}_p \ , \ \ p=1,...,m
		\end{equation*}
		we get a diagonal matrix of eigenvalues
		\begin{equation*}
			\Lambda \ = \
			\begin{pmatrix}
				\lambda_1\\
				&\ddots\\
				&&\lambda_m\\
			\end{pmatrix}
		\end{equation*}
		and the corresponding matrix {\bf R} containing the eigenvectors $\mathbf{x_p}$ in each column
		\begin{equation*}
			\mathbf{R}\ = \left( \mathbf{x_1} | \mathbf{x_2} | \dots | \mathbf{x_p} \right) 
		\end{equation*}
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Wave equation as Hyperbolic System}
		The Jacobian matrix \textbf{A} can now be expressed with the definitions  
		\begin{equation*}
			\begin{split}
				\textbf{A} &= \textbf{R}\Lambda\textbf{R}^{-1} \\
				\Lambda &= \textbf{R}^{-1}\textbf{A}\textbf{R} \ . 
			\end{split}
		\end{equation*}
		Applying these definitions to the wave equation
		\begin{equation*}
			\textbf{R}^{-1}\partial_t Q + \textbf{R}^{-1}\textbf{R}\Lambda\textbf{R}^{-1} \partial_x Q = 0 
		\end{equation*}
		and introducing the solution vector $\textbf{W}=\textbf{R}^{-1}\textbf{Q}$ results in
		\begin{equation*}
			\partial_t  \textbf{W} + \Lambda \partial_x \textbf{W} = 0 \ .
		\end{equation*}	
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Wave equation as Hyperbolic System}
		With $\lambda_{1,2}=\sqrt{\mu/\rho}= \pm c$, where c corresponds to the shear velocity, the eigenvectors can be obtained
		\begin{equation*}
			\textbf{r}_{1,2} \ = \
			\begin{pmatrix}
				\pm \rho c \\
				1
			\end{pmatrix}
		\end{equation*}
		interestingly enough containing as elements values of the {\it seismic impedance} $Z=\rho c$ that are relevant for the reflection behaviour of seismic waves. Thus the matrix \textbf{R} (and its inverse) are
		\begin{equation*}
			\textbf{R} \ = \ 
			\begin{pmatrix}
				Z & -Z \\
				1 & 1 
			\end{pmatrix}
			 \ , \ 
			\textbf{R}^{-1} \ = \ \frac{1}{2Z} \ 
			\begin{pmatrix}
				1 & Z \\
				-1 & Z 
			\end{pmatrix} \ \ .
		\end{equation*}		
	\end{frame}

% ---------------------------------------------------------

	\begin{frame}{The Wave equation as Hyperbolic System}
		The wave equation in the rotated eigensystem can be stated as
		\begin{equation*}
			\partial_t
			\begin{pmatrix}   
				w_1 \\
				w_2 
			\end{pmatrix}
			\ + \ 
			\begin{pmatrix}   
				-c & 0  \\
				0  & c
			\end{pmatrix}
			\partial_x
			\begin{pmatrix}   
				w_1 \\
				w_2 
			\end{pmatrix}
			\ = \ 0
		\end{equation*}
		with the simple general solution $w_{1,2} = w^{(0)}_{1,2}(x \pm ct)$, where the upper index 0 stands for the initial condition (i.e., waveform that is advected). The initial condition also fullfills  $\textbf{W}^{(0)}=\textbf{R}^{-1}\textbf{Q}^{(0)}$, we can therefore relate the so-called characteristic variables $w_{1,2}$ to the initial conditions of the physical variables as
		\begin{equation*}
			\begin{split}
				w_1(x,t)\ & = \ \frac{1}{2Z}(\sigma^{(0)}(x+ct) + Z   v^{(0)}(x+ct)) \\
				w_2(x,t)\ & = \ \frac{1}{2Z}(-\sigma^{(0)}(x-ct) + Z   v^{(0)}(x-ct))
			\end{split}
		\end{equation*} 
	\end{frame}		
		
% ---------------------------------------------------------

	\begin{frame}{The Wave equation as Hyperbolic System}
		To obtain the final analytical solution for velocity $v$ and stress $\sigma$ using $\textbf{Q}=\textbf{R}\textbf{W}$ as
		\begin{equation*}
			\begin{split}
				\sigma(x,t)\ =& \ \frac{1}{2}(\sigma^{(0)}(x+ct) +  \sigma^{(0)}(x-ct)) \\
             	& \ + \frac{Z}{2}(v^{(0)}(x+ct) -  v^{(0)}(x-ct)) \\
				v(x,t)    =& \ \frac{1}{2Z}(\sigma^{(0)}(x+ct) -  \sigma^{(0)}(x-ct)) \\
            	 & \ + \frac{1}{2}(v^{(0)}(x+ct) +  v^{(0)}(x-ct))  \ .
			\end{split}
		\end{equation*} 
		In compact form, the solution can be expressed as
		\begin{equation*} 
			\textbf{Q}(x,t) \ = \ \sum_{p=1}^m w_p(x,t) r_p 
		\end{equation*}
	\end{frame}
		
% ---------------------------------------------------------

\subsection{Solution to Scalar Advection Equation}

	\begin{frame}{Solution to Scalar Advection Equation}
		\small
		Denoting $q(x,t)$ for the unknown {\it scalar} solution and $a$  for the given (possibly space dependent) advection (wave) velocity to obtain the source-free wave equation
		\begin{equation*}
			 \partial_{t} q (x,t)  + a \  \partial_x q(x,t) = 0 	
		\end{equation*}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_dg_elements.png}
		\end{figure}
		The physical domain of each element is denoted by $D_k$ with left and right boundaries $x^k_l$ and $x^k_r$, respectively. Note the varying element sizes (h-adaptivity).
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Solution to Scalar Advection Equation}
	\small
		The k elements are \textbf{non-overlapping} which implies a higher number of degrees of freedom. The solution vector in 1D is
		\begin{equation*}
			q_{N_p} = \left(
				\begin{array}{cccc}
					q^1_1 & q^2_1 & \dots & q^n_1  \\
					q^1_2 & q^2_2 & \dots & q^n_2  \\
					\vdots & \vdots & \ddots & \vdots  \\
					q^1_{N_p} & q^2_{N_p} & \dots & q^n_{N_p} 
				\end{array}
			\right)
		\end{equation*}
		where $N_p$ is the number of points per element and $n$ the overall number of elements. The collocation points are Gauss-Lobatto-Legendre stored as
		\begin{equation*}
			x_{N_p} = \left(
			\begin{array}{cccc}
				x^1_1 & x^2_1 = x_{N_p}^1 & \dots & x^n_1  \\
				x^1_2 & x^2_2 & \dots & x^n_2  \\
				\vdots & \vdots & \ddots & \vdots  \\
				x^1_{N_p} = x_1^2 & x^2_{N_p} & \dots & x^n_{N_p} 
			\end{array}
			\right) \ .
		\end{equation*} 
	\end{frame}
		
% ---------------------------------------------------------

	\begin{frame}{Solution to Scalar Advection Equation}
		The global solution requires the adoption of flux concepts. Therefore a special sign is introduced that links the element-level to the global solution as
		\begin{equation*}
			q(x,t) \approx q_h(x,t) = \bigoplus_{k=1}^n q^k_h(x,t) 
		\end{equation*}
		where $\bigoplus$ indicates global assembly, $n$ is the overall number of elements and $h$ is the size of element $k$.
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
		To obtain the weak formulation of the scalar advection equation we multiply it 
		\begin{equation*}
			 \partial_{t} q (x,t)  + a \  \partial_x q(x,t)  = 0    
		\end{equation*}
		with a general test function $\phi_j(x)$, integrate over the k-th element domain $D_k$ to obtain
	\begin{equation*}
			 \int_{D_k} \partial_t q (x,t) \phi_j(x) dx + \int_{D_k} a \partial_x q (x,t) \phi_j(x) dx  = 0  
	\end{equation*}
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
	And partially integrate replacing the term containing the space derivative by
	\begin{equation*}
		\begin{split}
			 \int_{D_k} a \partial_x q (x,t) \phi_j(x) dx \ =& \ [a q(x,t) \phi_j(x) ]_{x_l}^{x_r}  \\
			 -&\int_{D_k} a q (x,t) \partial_x \phi_j(x) dx \ 
		\end{split}
	\end{equation*}
where $x_r$ and $x_l$ are the right and left boundaries of element $k$, and we assume constant velocity $a$ inside the element.  
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
	\begin{minipage}{0.3\textwidth}
		\begin{figure}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_dg_gauss.png}
		\end{figure}
	\end{minipage}
	\mbox{}		
	\hfill		
	\begin{minipage}{0.65\textwidth}
			What does the rule of partial integration look like with more than one dimension? The definition is known as
		\begin{equation*}
			\begin{split}
				\int_{\Omega} \partial_{x_i} u v d\Omega =& \int_{\Gamma} u v n_i d\Gamma  \\
				-& \int_{\Omega} u \partial_{x_i} v d\Omega 
			\end{split}
		\end{equation*}
		where $u,v$ are arbitrary space-dependent functions, $\Omega$ denotes the entire volume, $\Gamma$ its boundary, and $n_i$ is a vector normal to the boundary.
	\end{minipage}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
		Putting this into the advection equation we obtain for element k
		\begin{equation*}
			\begin{split}
				\int_{D_k} \partial_t q (x,t) \phi_j(x) dx \ -& \ \int_{D_k} a q (x,t) \partial_x \phi_j(x) dx \\
				=& \ - [a q(x,t) \phi_j(x) ]_{x_l}^{x_r} \ 
			\end{split}
		\end{equation*}
		Replacing the unknown field $q(x,t)$ by a finite polynomial representation in terms of a weighted sum over Lagrange polynomials inside each element k of order N denoted as $\ell_i (x), i=1,\dots,N+1$ defined in the interval $x\in [-1,1]$. 
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
		With the exact interpolation property of the Lagrange polynomials at the Gauss-Legendre-Lobatto collocation points $x_i$, here illustrated for an arbitrary function $y_i = y(x_i)$ 
		\begin{equation*}
			y_i=\sum_{j=1}^{N_p} y_j \ell_j(x_i)=\sum ^{N_p}_{j=1}y_j \delta_{ij}
		\end{equation*}
		we obtain for element $k$      
		\begin{equation*}
			 q(x,t) |_{x=x_i} \ = \sum_{i=1}^{N_p} q_i(t) \ \ell_i(x)
		\end{equation*}
		where $N_p$ indicates that in each element the polynomial order may vary, $x=x_i$ are the collocation points.
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
		In addition we also replace the test function by Lagrange polynomials. Combining this we yield
		\begin{equation*}
			\int_{D_k} \partial_t \sum_{i=1}^{N_p} q_i (t) \ell_i(x) \ell_j(x) dx - \int_{D_k} a \sum_{i=1}^{N_p} q_i (t) \ell_i (x) \partial_x \ell_j(x) dx
		\end{equation*}
		and after re-ordering (omitting space dependencies of polynomials)
		\begin{equation*}
			\sum_{i=1}^{N_p} \left[ \left[ \int_{D_k}  \ell_i(x) \ell_j(x) dx \right]  \partial_t q_i(t)  - \left[ \int_{D_k} a \ell_i(x) \partial_x \ell_j(x) dx \right] q_i (t)  \right]
		\end{equation*}
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Weak formulation}
	\small
		Recognizing the elemental mass $M_{ij}$ and stiffness $K_{ij}$ matrices in this equation. Assuming implicit matrix-vector operations we obtain
		\begin{equation*}
			M \ \partial_t q(t) \ - \ K^T  \ q(t)
		\end{equation*}
		where the matrices are given by
		\begin{equation*}
			\begin{split}
				M_{ij} \ &= \ \int_{D_k}  \ell_i(x) \ell_j(x) dx \\
				K_{ij} \ &=  \ \int_{D_k} a  \ell_i(x) \partial_x \ell_j(x) dx
			\end{split}
		\end{equation*}
		assuming constant velocity $a$ inside element $k$. Note that the lower index $D_k$ indicates that we are still in physical space and we have to map to the local coordinate system.
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Elemental Mass and Stiffness Matrices}
		The matrices (nodal or modal approach) for arbitrary test functions $\phi_i(x)$ are defined by
		\begin{equation*}
			\begin{split}
				M_{ij}\ &= \ \int_{D_k} \phi_i(x) \phi_j(x) \ dx \\
				K_{ij}\ &= \ a \ \int_{D_k} \phi_i(x) \partial _x  \phi_j(x) \ dx \\
			\end{split}
		\end{equation*}
		containing integrals over (derivatives of) the test functions $\phi(x)$. Replacing the integrals by a weighted sum over the function values $f(x_i)$ at carefully chosen points $x_i$ inside the elements
		\begin{equation*}
			\int_\Omega f(x) \ dx \approx \sum_{i=1}^N w_i \ f(x_i) 
		\end{equation*}
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{Elemental Mass and Stiffness Matrices}
		We need to map our physical coordinates into an element-based system.  In 1D this is quite simple using as local variable $\xi$ and transforming via
		\begin{equation*}
			x^k(\xi) = x_l^k + \frac{(1+\xi)}{2} dx \ \ \xi \in [-1,1]
		\end{equation*}
		where $x_l^k$ and $x_r^k$ are the left and right physical boundaries of element $k$. In general the mapping of the differential used to evaluate integrals is called the Jacobian defined for element $k$ as   
		\begin{equation*}
			J^k=  \frac{dx}{d \xi} \ , \ \ \ J^k=\frac{x^k_r-x^k_l}{1- (-1)} = \frac{h^k}{2}  
		\end{equation*}
		where $h^k$ is the size of element $k$. 
	\end{frame}	
						
% ---------------------------------------------------------

	\begin{frame}{Elemental Mass and Stiffness Matrices}
		For the elemental matrices we obtain for general test functions
		\begin{equation*}
			\begin{split}
				M_{ij}\ &= \ \int_{D_k} \phi_i(\xi) \phi_j(\xi) J^k \ d\xi \\
				K_{ij}\ &= \ a \ \int_{D_k} \phi_i(\xi) \partial_\xi  ((J^k)^{-1}) J^k \phi_j(\xi) \ d\xi \\
				&= \  a \int_{D_k} \phi_i(\xi) \partial_\xi  \phi_j(\xi) d\xi 
			\end{split}
		\end{equation*}
	\end{frame}
							
% ---------------------------------------------------------

	\begin{frame}{Elemental Mass and Stiffness Matrices}
		Finally, we can replace the test function with the Lagrange polynomials of order $N_p$ leading to the definition of the mass and stiffness matrices
		\begin{equation*}
			\begin{split}
				M_{ij}=\int_{-1}^1 \ell_i(\xi) \ell_j(\xi) \ J^k \ d\xi =& \sum_{m=1}^{N_p} w_m \ \ell_i(x_m)  \ell_j(x_m)\ J^k \\
				=&\sum_{m=1}^{N_p} w_m \delta_{im}\ \delta_{jm} \ J^k \\ 
				=& \begin{cases} w_i \ J^k \ \ \text{ if } i=j \\ 0 \ \ \ \ \ \ \ \text{   if } i \neq j\end{cases}
			\end{split}	
		\end{equation*}
	\end{frame}		
							
% ---------------------------------------------------------

	\begin{frame}{Elemental Mass and Stiffness Matrices}		
		\begin{equation*}
			\begin{split}
				K_{ij}= \int_{-1}^1 \ell_i(\xi) \partial _x  \ell_j(\xi) \ d\xi=& \sum_{m=1}^{N_p} a \ w_m \ \ell_i(x_m) \partial_x \ell_j(x_m) \\
				=& \sum_{m=1}^{N_p} a \ w_m \delta_{im}\partial_x \ell_j(x_m) \\
				=& a \ w_i \partial_x \ell_j(x_i) \ . 
			\end{split}
		\end{equation*}
	\end{frame}
								
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		The key task in the discontinuous Galerkin scheme concerns the question what values to allocate to the points at the element boundaries. This involves the use of flux schemes.
		\begin{figure}
			\includegraphics[width=0.45\textwidth]{\figpath/fig_dg_flux.png}
		\end{figure}
	\end{frame}	
									
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		Starting with the r.h.s. of this equation
		\begin{equation*}
			\begin{split}
				\int_{D_k} \partial_t q (x,t) \phi_j(x) dx \ -& \ \int_{D_k} a q (x,t) \partial_x \phi_j(x) dx \\
				=& \ - [a q(x,t) \phi_j(x) ]_{x_l}^{x_r} \ 
			\end{split}
		\end{equation*}
		The originial form holds also for higher dimensions
		\begin{equation*}
		 \int_{\partial D_k} a  \  q(x,t) \phi_j(x) \ \mathbf{n} \ dx
		\end{equation*}
		where $\mathbf{n} = \pm 1$ denotes the vector normal to the boundary, in 1D taking the values $n=-1$ and $n=1$ at the left and right boundaries.
	\end{frame}
										
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
	\small
		Replacing the space-dependent part of $q(x,t)$ by a sum over Lagrange polynomials to obtain   
		\begin{equation*}
			\sum_{i=1}^{N_p} \int_{\partial D_k} \ell_j(x) \ \ell_i(x) \ a \ q_i(t) \ \mathbf{n} \ dx
		\end{equation*}
		and also replacing the arbitrary test function with the same function. The orthogonality of the Lagrange polynomials and the fact that we are integrating over the surface $\partial D_k$ leads to
		\begin{equation*}
			\sum_{i=1}^{N_p}(\ell_i(x_r^k) \ \ell_j(x_r^k) (a \ q(x_r^k))^* -\ell_i(x_l^k) \ \ell_j(x_l^k)(a \ q(x_l^k))^*) 
		\end{equation*}
		\begin{equation*}
			= \ell_j(x_r^k) (a \ q(x_r^k))^* - \ell_j(x_l^k)(a \ q(x_l^k))^*
		\end{equation*}
		where we introduced the {\it starred} terms $(a \ q(x_{r,l}^k))^*$ that are the values at the element boundaries.
	\end{frame}
											
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		Therefore, we introduce the general flux vector with elements $F_j(a,q^k(x,t)), j=1,..,N_p$ as
		\begin{equation*}
			F_j(a,q^k(x,t)) \ = \ [\ell_j(x) \ (a \ q(x,t))^*]_{x^{k}_l}^{x^{k}_r} 
		\end{equation*}
		In vector form
		\begin{equation*}
\textbf{F} = \left(
\begin{array}{c}
F_1 \\
0 \\
\vdots \\
0\\
F_{N_p} 
\end{array}
\right)
\ = \ 
 \left(
\begin{array}{c}
-(a \ q(x,t))^*(x^k_l) \\
0 \\
\vdots \\
0 \\
(a \ q(x,t))^*(x^k_r)
\end{array}
\right) \ . 
\end{equation*}
	\end{frame}
												
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_dg_bc.png}
		\end{figure}
		At the left boundary the outside value in element $k-1$ is denoted with "`+"' and the value inside element $k$ with "`-"'. The starred value $q*(x_l^k)$ is the flux that needs to be defined as a function of the boundary values of the adjacent elements.
	\end{frame}	
													
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		Taking the average of the values from both sides of the boundaries - also called the \textit{central flux} - can be expressed as
		\begin{equation*}
			\begin{split}
				F^c_1 \ &=  \   \frac{1}{2} a^k (q(x^{k-1}_r,t) + q(x^k_l,t)) \\
				F^c_{N_p} \ &=  \   \frac{1}{2} a^k (q(x^{k}_r,t) + q(x^{k+1}_l,t))
			\end{split}
		\end{equation*}		
		The simplemost, stable choice is the so called \textit{upwind flux}
		\begin{equation*}
			\begin{split}
				F^{up}_1 \ &= \ \begin{cases} a^k \ q(x^k_l) \ \ \ \ \text{if  } \ \ a^k \leq 0\\ a^{k-1} \ q(x^{k-1}_r) \ \text{if  } \ \ a^{k-1} > 0\end{cases} \\
				F^{up}_{N_p} \ &= \ \begin{cases} a^{k+1} \ q(x^{k+1}_l) \ \ \ \ \text{if  } \ \ a^{k+1} \leq 0\\ a^k \ q(x^k_r) \ \text{if  } \ \ a^k > 0\end{cases}
			\end{split}
		\end{equation*}
	\end{frame}
														
% ---------------------------------------------------------

	\begin{frame}{The Flux Scheme}
		Centered and upwind fluxes can be formulated in a compact way. This formulation reads
		\begin{equation*}
			\begin{split}
				F_1 &= -a \frac{1}{2}(q(x^k_l) + q(x^{k-1}_r))  -  \frac{|a|}{2} (1-\alpha) (q(x^{k-1}_r)-q(x^k_l)) \\
				F_{N_p} &=    a  \frac{1}{2}(q(x^k_r) + q(x^{k+1}_l))  +  \frac{|a|}{2} (1-\alpha) (q(x^k_r) - q(x^{k+1}_l))
			\end{split}
		\end{equation*}
		where $\alpha=0$ corresponds to the upwind flux and $\alpha=1$ to the centered flux scheme. 
	\end{frame}
															
% ---------------------------------------------------------

\section{The Discontinuous Galerkin Method Put To Action}

																
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		In matrix notation we yielded for one element
		\begin{equation*}
			 M \partial_{t} q(t) \ - \  K^T q(t) \  = \  -F(a,q(t)) 
		\end{equation*}
		requiring an extrapolation scheme of the form
		\begin{equation*}
			\partial_{t} q(t) \  = \ M^{-1} (K^T q(t)  \ -F(a,q(t)))  
		\end{equation*}
		where $F(a,q(t))$ is the flux vector as defined above. 
		We seek to extrapolate the system from some initial conditions and obtain using the simple Euler method for each element
		\begin{block}{}
		\begin{equation*}
			q(t_{n+1})  \ \approx \ q(t_{n}) \ + dt \left[ M^{-1}( K^T q(t_n) \ - \ F(a,q(t)) \right]
		\end{equation*}
		\end{block}
		where for the flux scheme F() we use the upwind approach.
	\end{frame}
																	
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		Employing a high-order extrapolation procedure known as \textit{predictor-corrector} method (or Heun's method, or two-stage Runge-Kutta method). At time step $t_i$ using time increment dt
		\begin{equation*}
			\begin{split}
				k_1&=f(t_i,y_i)\\
				k_2&=f(t_i+dt,y_i+dt k_1)\\
				y_{i+1}&=y_i+\frac{1}{2}dt (k_1+k_2) \ 
			\end{split}
		\end{equation*}
	\end{frame}
																		
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_dg_matrix.png}
		\end{figure}
		The matrix-vector form of the discontinuous Galerkin method. The system of questions at an elemental level is illustrated by plotting the absolute matrix/vector values. N = 3, $N_p$ = N + 1 = 4
	\end{frame}	
																	
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		
		\begin{minipage}{0.45\textwidth}
		\scriptsize
		\begin{mdframed}[style=MyFrame]
{\sffamily
{\scriptsize
\% Initialize vectors, matrices \\
Minv = zeros(N+1,N+1,ne); \\
K = zeros(N+1,N+1,ne);\\
q = zeros(N+1,ne);\\
(...)\\
for k=1:ne,\\
\indent for i=1:N+1,\\
\indent \indent     Minv(i,i,k)=1./(w(i)*J(k); \\
\indent end \\
end \\
(...)\\
for k=1:ne, \\
\indent for i=1:N+1, \\
\indent  \indent for j=1:N+1, \\  
\indent \indent \indent    K(i,j,k)=a(k)*w(i)*l1d(j,i); \\
\indent \indent end \\
\indent end \\
end
}}
\end{mdframed}
		\end{minipage}
	\mbox{}		
	\hfill			
		\begin{minipage}{0.5\textwidth}
			\small
			The most important initialization step is the calculation of the elemental matrices, mass $M$ and stiffness $K$. The following code part illustrates a possible implementation looping through all elements $ne$ and calculating these matrices as a function of the Jacobian $J(k)=h^k/2$ where $h^k$ is the size  of element $k$:   
		\end{minipage}
	\end{frame}
																		
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		
		{\scriptsize
		\begin{mdframed}[style=MyFrame]
{\sffamily
{\scriptsize
\% Flux calculation \\
F = zeros(N+1,ne);\\
(...)\\
for k=1:ne,\\
F(1,k)= -a*(q(1,k)+q(N+1,k-1))/2 - abs(a)*(1-alpha)/2*(q(N+1,k-1)-q(1,k)); \\
F(N+1,k)= a*(q(N+1,k)+q(1,k+1))/2 + abs(a)*(1-alpha)/2*(q(N+1,k)-q(1,k+1)); \\
end
}}
\end{mdframed}
}
\vspace{1cm}
		The flux vector has to be calculated new for each time step (or intermediate step when using high-order extrapolation schemes) as it depends on the current values of the wavefield q. $alpha$ can be used to change the flux scheme from central ($alpha=1$) to upwind ($alpha=0$).
	\end{frame}
																			
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
		The element-wise system of equations can be extrapolated by the Euler scheme as
		\vspace{0.3cm}
		{\scriptsize
		\begin{mdframed}[style=MyFrame]
{\sffamily
{\scriptsize
\% Time extrapolation \\
for it = 1:nt, \\
(...)\\
\% Initialize flux vectors F
(...) \\
\indent for k=1:ne, \\
\indent \indent q(:,k) = q(:,k) + dt .* (Minv(:,k) * (K(:,k)'*q(:,k)-F(:,k))); \\
\indent end
(...)\\
end
}}
\end{mdframed}
}
\vspace{0.3cm}
where $nt$ is the overall number of time steps, $dt$ is the global time increment, and $ne$ is the number of elements. 
Note that can directly update the solution vector $q$ without intermediate storage at different time level(s). 
	\end{frame}
																			
% ---------------------------------------------------------

	\begin{frame}{The Discontinuous Galerkin Method Put To Action}
			A high-order extrapolation scheme like the predictor-corrector method can be implemented as\\
			\vspace{0.2cm}
			\begin{minipage}{0.45\textwidth}
				\scriptsize 		
			\begin{mdframed}[style=MyFrame]
{\sffamily
{\scriptsize
\% Time extrapolation \\
for it = 1:nt, \\
(...)\\
\% Predictor corrector scheme \\
\% Initialize flux vectors F for all k \\
(...) \\
\% First step (predictor) \\
\indent for i=1:k,\\
\indent \indent k1(:,k)= Minv(:,k) * (K(:,k)'*q(:,k)-F(:,k));\\
\indent end\\
(...) \\
}}
\end{mdframed}
			\end{minipage}
	\mbox{}		
	\hfill					
			\begin{minipage}{0.48\textwidth}
			 \scriptsize
			\begin{mdframed}[style=MyFrame]
{\sffamily
{\scriptsize
\% Initialize flux vectors F for q+dt*k1\\
(...) \\
\% Second step (corrector)\\
\indent for i=1:k,\\
\indent \indent k2(:,k)= Minv(:,k) * (K(:,k)'*(q(:,k)+dt*k1(:,k))-F(:,k));\\
\indent end\\
(...)\\
\% Update\\
\indent for i=1:k,\\
\indent \indent q(:,k)= q(:,k) + dt (k1(:,k) + k2(:,k));\\
\indent end\\
(...)\\
end
}}		
\end{mdframed}	
			\end{minipage}
	\end{frame}	
																				
% ---------------------------------------------------------

	\begin{frame}{Example}
	\footnotesize
		\begin{table}
\caption{\footnotesize Simulation parameters for 1D discontinuous Galerkin advection }
\begin{tabular}[t]{lll}
Parameter & Value & Meaning \\
\hline
$ne$ & 200 & Elements\\
$N $ & 3 & Order \\
$ a$ & 2500 $m/s$ & Velocity\\
$x_{max}$ & 10000 m & x-domain\\
$dx_{min}$ & 13.82 m & Increment\\
dt & 4.4e-4 s & Time step\\
eps & 0.08 & Courant \\
$\sigma$ & 300 m & Gauss width\\
$x_0$ & 1000 m & Source x\\
$t_{max}$ & 3 s & Duration
\end{tabular}
		\end{table}
		Initiating the simulation with a spatial initial condition using a Gaussian function $e^{-1/\sigma^2(x-x_0)^2}$. A source term can be added to the system of equations in a straight forward way.
	\end{frame}	
																				
% ---------------------------------------------------------

	\begin{frame}{Example}
		\begin{figure}
			\includegraphics[width=0.65\textwidth]{\figpath/fig_dg_comp.png}
		\end{figure}
	\end{frame}	
																				
% ---------------------------------------------------------

	\begin{frame}{Summary}
	\scriptsize
		\begin{itemize}
			\item The discontinuous Galerkin method is a finite-element type method. The main difference is that the solution fields are not continuous at the element boundaries.
\item The elemental mass and stiffness matrices are formulated very similar to the classic finite-element schemes. However, they are never assemble to a global system of equations. Therefore no large system matrix needs to be inverted.
\item Elements are linked by a flux scheme, similar to the finite volume method. This scheme leads to an entirely local algorithm in the sense that all calculations are carried out at an elemental level. Communication happens only to direct neighbours.
\item The discontinuous Galerkin scheme can easily to higher orders, keeping the local nature of the solution scheme. This leads to high efficiency when parallelizing the scheme.
		\end{itemize}
	\end{frame}	
																					
% ---------------------------------------------------------

	\begin{frame}{Summary}
	\scriptsize
		\begin{itemize}
\item The solution fields can be extended using nodal and modal approaches. Modal approaches are preferred for tetrahedral schemes. Nodal solutions are preferred for hexahedral meshes.
\item The discontinuous behaviour at the element boundaries and the associated discretization of the element boundaries increases the number of degrees of freedom compared to other methods.
\item The flexibility with polynomial order, element size, local time stepping leads to a formidable problem when parallelizing a discontinuous Galerkin method: load balancing the computational task is not easy! Solution to this problem requires tight cooperation with computational scientists.
\item In seismology, the discontinuous Galerkin method is useful for problem with highly complex geometries (by using tetrahedral meshes) and for problem with non-linear internal boundary conditions (e.g., dynamic ruputure problems).
		\end{itemize}
	\end{frame}	
																					
% ---------------------------------------------------------

	\begin{frame}{Comprehensive Questions}
	\scriptsize
		\begin{enumerate}
\item
What were the key  points that led to the development of the discontinuous Galerkin method in seismology? Discuss the pros and cons of the method compared to finite-element type methods and the finite-difference method. 
\item
Explain qualitatively the difference between nodal and modal approaches. 
\item
Explain why the discontinuous Galerkin method lends itself to parallel implementation on supercomputer hardware. 
\item 
What is $p-$ and $h-$adaptivity? Why is it straight forward to have this adaptivity with the discontinuous Galerkin method and not with others? Give examples  in seismology where this adaptivity can be exploited and why. 
\item
What is local time-stepping? For what classes of Earth models and/or problems in seismology might it be useful?
\item
What is the problem that arises on computers when using algorithms with h-/p-adaptivity and local time stepping? 
\item
What is meant by conservative and non-conservative properties in the context of advection problems?
\end{enumerate}
	\end{frame}	
																						
% ---------------------------------------------------------

	\begin{frame}{Theoretical Problems}
		\scriptsize
		\begin{enumerate}
\setcounter{enumi}{7}
\item
Show that the advection problem $\partial_t q + a \partial_x q = 0$ has a hyperbolic form. 
\item
The coupled 1D wave equation for longitudinal velocity $v$ and pressure $p$ can  be formulated with compressibility $K$ and density $\rho$ as  
\begin{equation*}
\begin{split}
\partial_t p + K \partial_x v &= 0 \\
\partial_t v + \frac{1}{\rho} \partial_x p &= 0 \ .
\end{split}
\end{equation*}
Formulate the Jacobian matrix of the coupled system of equations and calculate its eigenvalues. 
\item
Show that the rule of partial integration corresponds to Gauss' theorem in higher dimensions (assuming on of the functions under the integral to be unity). Explain the relevance of this for the discontinuous Galerkin method. 
\item
Reformulate the nodal discontinuous Galerkin solution allowing the velocity to be variable inside the element (Note: Compare with the spectral-element formulation). 
\item
Show that setting $\alpha=0$ leads to the upwind flux scheme. 
\item
Search in the literature for the {\it classical} 4-term Runge-Kutta method. Formulate a pseudo-code for the scalar advection problem for this extrapolation scheme. 
\end{enumerate}
	\end{frame}
																							
% ---------------------------------------------------------

	\begin{frame}{Programming Exercises}
		\scriptsize
		\begin{enumerate}
\setcounter{enumi}{13}
\item
Using the sample code $dg1d$ find numerically the Courant limit (for a fixed time extrapolation to $t_{max}$), keeping all other parameters constant, increasing the global spatial order $N$ of the scheme.
\item
Formulate an upwind finite-difference scheme for the scalar advection problem and write a computer programm. Discuss the diffusive behaviour. 
\item
Modify the sample code $dg1d$ such that each element can have its own polynomial order ($p$-adaptivity) and size ($h-$adaptivity). (Suggestion: Initialize the shape of the solution and other matrices using the maximum number of degrees of freedom $N^p_{max}$).     
\item
Extend the sample code $dg1d$ for the scalar advection problem to the 4-term Runge-Kutta method. Compare the accuracy of the method with the lower order extrapolation scheme as a function of spatial order $N$ inside the elements. 
\end{enumerate}
	\end{frame}
																							
% ---------------------------------------------------------

	\begin{frame}{Programming Exercises}
		\scriptsize
		\begin{enumerate}
\setcounter{enumi}{17}	
\item
Formulate the analytical solution to the advection problem and plot it along with the numerical solution in $dg1d$ each time you visualize it during extrapolation. Formulate an error between analytical and numerical result. Analyze the solution error as a function of propagation distance for the Euler scheme and the predictor-corrector scheme. 
\item
Explore the $p-$ and $h-$ adaptivity of the discontinuous Galerkin method $\ $ in the following way. Using an appropriate Gaussian function defined on the entire physical domain decrease the element size by a factor of 5 towards the center of the domain. Find an appropriate variation of the order inside the elements to obtain a reasonable computational scheme (in the sense that the grid point distance does not vary too much). Hint: Use high order schemes at the edges of the physical domain and low(est) order schemes at the centre of the domain.   
\end{enumerate}
	\end{frame}
	
	
	





% ---------------------------------------------------------

\end{document}
