

% beamer version of chapter 4 "The Finite Volume Method" in "Computational Seismology,  A Practical Introduction" (H. Igel)

\documentclass[aspectratio=169]{beamer}

\mode<presentation>
{
  %\usetheme{metropolis}
	\usetheme{metropolis}
  %\usecolortheme{beaver}
	\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}
\usepackage{mathrsfs}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%Implement Matlab and Python Codes%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\usepackage{listings}
\lstloadlanguages{bash,  Python,  Matlab}
\lstset{
alsolanguage=bash, %Setzen der Codesprache
alsolanguage=Python,
alsolanguage=Matlab,
basicstyle=\tiny \sffamily , % alle listings winzig drucken
keywordstyle=\color{blue} \bfseries , % Schlüsselwörter fett und schwarz drucken
identifierstyle=\color{black},
stringstyle=\color{orange}, % Strings werden orange gedruckt
%otherkeywords={*,+,-}=\color{green}, % andere Schlüsselwörter, die gekennzeichnet werden sollen
commentstyle=\color{purple!40!black}, % Kommentare blau drucken
%stringstyle=\ttfamily, % Strings im Code Schreibmaschinenähnlich - setzt sich etwas vom Code ab
backgroundcolor=\color{white!95!black},
showstringspaces=false, % Strings im Code ohne Kenntlichmachung von Leerzeichen
captionpos=b,   %Caption erscheint unterhalb des Codes
frame=single, %Code wird umrandet
numbers=none, % Nummerierung der Linien links vom Code
numberstyle=\footnotesize, % Art der Nummerierung, hier z. B.: winzige Schrift
%stepnumber=5, % Intervall der Nummern, hier z.B. 5
              % falls nicht gesetzt: Alle Nummern werden ausgegeben.
numbersep=9 pt % Abstand vom Code, hier z.B. 5 pt
} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\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_fv}
\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Finite Volume 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.45\textwidth}
			\begin{figure}
				\includegraphics[width=0.95\textwidth]{\figpath/fig_fv_mesh1.png}
			\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill		
		\begin{minipage}{0.5\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item Robust
					\item Simple concept
					\item Irregular grids
					\item Explicit method
					\item Find solution for strongly heterogeneous paramters and discontinuities
				\end{itemize}	
			\end{block}	
		\end{minipage}
	\end{frame}

% ---------------------------------------------------------
\subsection{History}
	
	\begin{frame}{History}
           \begin{itemize}
           	\item First applications in plasma physics (Hermeline, 1993) and computational fluids (Versteeg and Malalasekera, 1995)
           	\item Discrete version of the divergence theorem was used for seismic wave propagation (Dormy and Tarantola, 1995)
           	\item Comparing and quantifying the accuracy of wave propagation of the finite volume method was done by Kaeser et al. 2001
           	\item Leveque, 2001 presented a natural consequence of conservation laws
           	\item Dumbser et al., 2007a presented the arbitrary high-order scheme (ADER) to the finite-volume method
			\end{itemize}	
	\end{frame}
	
% ---------------------------------------------------------	
	
	\begin{frame}{History}
			\begin{figure}
				\includegraphics[width=0.55\textwidth]{\figpath/fig_fv_tets.png}
			\end{figure}		
	\end{frame}

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

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

\subsection{Finite Volume in a Nutshell}	
	
	\begin{frame}{Finite Volume in a Nutshell}
		\begin{minipage}{0.6\textwidth}
			In its basic form it takes an entirely local viewpoint in the sense that the solution field q(x,t) is tracked inside a cell. The field is approximated by an average quantity $Q_i^n$ inside cell C as:
			\begin{equation*}
				Q^n_i \ = \ \frac{1}{dx} \int_C q(x,t) dx
			\end{equation*}
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.35\textwidth}
			\begin{figure}
				\includegraphics[width=0.65\textwidth]{\figpath/fig_fv_nutshell.png}
			\end{figure}
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Finite Volume in a Nutshell}
		Tracking the change of the values with time inside each cell
		\begin{equation*}
			\frac{Q_i^{n+1} - Q_i^n}{dt} \ = \ \frac{F^n_{i-1/2} -  F^n_{i+1/2}}{dx} 
		\end{equation*}
		Assuming the flux depends only on the adjacent $Q_i^n$ values (for left boundary)
		\begin{equation*}
			F^n_{i-1/2} \ = \ f(Q_{i-1}^{n},Q_i^n) 
		\end{equation*}
		The requirement of conservation for a transport problem leads to the advection equation of the form
		\begin{equation*}
			\partial_t q(x,t) \ + \ a \ \partial_x q(x,t)  \ = \ 0
		\end{equation*}
		where a is the transport velocity
	\end{frame}

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

	\begin{frame}{Finite Volume in a Nutshell}
		For the constant-coefficient advection problem the flux terms obviously are
		\begin{equation*}
			\begin{split}
				F^n_{i-1/2} \ &= \ a Q^n_{i-1} \\
				F^n_{i+1/2} \ &= \ a Q^n_{i}
			\end{split}	
		\end{equation*}
		where a is the advection velocity. With these definitions  we obtain a fully discrete extrapolation scheme as 
		\begin{equation*}
			Q_i^{n+1} \ = \ Q_i^n + a \frac{dt}{dx} (Q_{i-1}^{n}-Q_i^n) 
		\end{equation*}
		A considerable better solution is the  Lax-Wendroff scheme given as
		\begin{equation*}
			Q_i^{n+1} \ = \ Q_i^n - \frac{a dt}{2 dx} (Q_{i+1}^{n}-Q_{i-1}^{n}) + \frac{1}{2}(\frac{a dt}{dx})^2 (Q_{i-1}^{n}-2 Q_{i}^{n}+ Q_{i+1}^{n}) 
		\end{equation*}
		which is second-order accurate and much less dispersive. 
	\end{frame}

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

\section{Ingredients}

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

\subsection{The Finite-Volume Method via Conservation Laws}
	\begin{frame}{The Finite-Volume Method via Conservation Laws}
		To describe what is happening we put ourselve into a finite volume cell that we denote as $\mathscr{C}$ and denote the boundaries as $x\in{x_1, x_2}$. We further assume a positive advection speed $a$. 
		\begin{figure}
			\includegraphics[width=0.4\textwidth]{\figpath/fig_fv_volume.png}
		\end{figure}
	\end{frame}

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

	\begin{frame}{The Finite-Volume Method via Conservation Laws}
		The total mass of any quantity inside the cell is 
		\begin{equation*}
			\int_{x_1}^{x_2} q(x,t) dx
		\end{equation*}
		and a change in time can only be due to fluxes across the left and/or right cell boundaries. Thus
		\begin{equation*}
			\partial_t \int_{x_1}^{x_2} q(x,t) dx \ = \ F_1(t) - F_2(t) 
		\end{equation*}
		where $F_i(t)$ are rates (e.g., in $g/s$) at which the quantity flows through the boundaries.
	\end{frame}

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

	\begin{frame}{The Finite-Volume Method via Conservation Laws}
		If we assume advection with a constant transport velocity $a$ this flux is given as a function of the values of $q(x,t)$ as
		\begin{equation*}
			F \rightarrow f(q(x,t))  \ = \ a q(x,t)
		\end{equation*}
		in other words
		\begin{equation*}
			\partial_t \int_{x_1}^{x_2} q(x,t) dx \ = \  f(q(x_1,t)) -  f(q(x_2,t))
		\end{equation*}
		This is called the integral form of a hyperbolic conservation law. 
	\end{frame}

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

	\begin{frame}{The Finite-Volume Method via Conservation Laws}
		Using the definition of integration and antiderivates to obtain 
		\begin{equation*}
			\begin{split}
				\partial_t \int_{x_1}^{x_2} q(x,t) dx \ = \  - \int_{x_1}^{x_2} \partial_x f(q(x,t)) dx  &\\
				\int_{x_1}^{x_2} \left[  \partial_t q(x,t) + \partial_x f(q(x,t)) \right] dx \ = 0 
			\end{split}
		\end{equation*}
		which leads to the well-known partial-differential equation of linear advection
		\begin{equation*}
		 	\partial_t q(x,t) + \partial_x f(q(x,t)) \ = \ 0 
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		\begin{minipage}{0.55\textwidth}
		Instead of working on the field q(x,t) itself we approximate the integral of q(x,t) over the cell $\mathscr{C}$ by
		\begin{equation*}
			Q^n_i \ \approx \ \frac{1}{dx} \int_\mathscr{C} q(x,t^n) dx 
		\end{equation*}
		This is the average value of $q(x,t)$ inside the cell.		
		\end{minipage}
    \mbox{ }
    \hfill	
		\begin{minipage}{0.4\textwidth}
		\begin{figure}
			\includegraphics[width=0.9\textwidth]{\figpath/fig_fv_flux.png}
		\end{figure}		
		\end{minipage}
	\end{frame}
	
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		In order to find an extrapolation scheme to approximate the future state of our finite-volume cells we integrate the integral form of a hyperbolic conservation law.
		\begin{equation*}
			\begin{split}
 				& \int_C q(x,t^{n+1}) dx  \ - \ \int_C q(x,t^{n}) dx \\
 				=& \int_{t_n}^{t_{n+1}} f(q(x_L, t) dt - \int_{t_n}^{t_{n+1}} f(q(x_R, t) dt
			\end{split}
		\end{equation*}
		where we used the definitions $\mathscr{C}\rightarrow [x_1, x_2] = [x_L, x_R]$, rearranged terms, and divided by $dx$ in order to recover the average cell values. \\
		\textbf{This equation is exact! }
	\end{frame}	
	
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		Using the following terms for the fluxes at the boundaries  
		\begin{equation*}
			F^n_{L,R} \ = \ \int_t^{t_{n+1}} f(q(x_{L,R}, t)) dt
		\end{equation*}
		we obtain a time-discrete scheme for the average values of our solution  field q(x,t)
		\begin{equation*}
		Q_i^{n+1} \ = \  Q_i^n \ - \ \frac{dt}{dx} (F^n_{R} - F^n_{L})
		\end{equation*}
		where the upper index $n$ denotes time level $t_n=n*dt$ and the lower index $i$ denotes cell $\mathscr{C}_i$ of size $dx$. 
	\end{frame}	
		
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_fv_upwind.png}
\caption{For the linear advection problem we can analytically predict where the tracer will be located after time $dt$. The value of $q(x_i,t^{n+1})$ will be exactly the same as  $q(x_i-a dt,t^n)$. We can use this information to predict the new cell average $Q_i^{n+1}$.}
		\end{figure}
	\end{frame}
			
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		We thus seek to approximate the next cell update $Q_i^{n+1}$ knowing that
		\begin{equation*}
			Q_i^{n+1} \approx  q(x_i,t^{n+1}) = q(x_i-a dt,t^{n}) 
		\end{equation*}
		The new cell average analytically by adding the appropriate mass flowing via the left boundary by interpolation 
		\begin{equation*}
			\begin{split}
				Q_i^{n+1} \ & = \ Q_{i-1}^n + \frac{dx - a dt}{dx} (Q_{i}^n-Q_{i-1}^n) \\
				Q_i^{n+1} \ & = \ Q_{i}^n (1 - \frac{a dt}{dx}) + Q_{i-1}^n \frac{a dt}{dx} \ . 
			\end{split}
		\end{equation*}
		After re-arranging we finally obtain a fully discrete scheme 
		\begin{equation*}
			Q_i^{n+1} \ = \ Q_i^{n} - \frac{a dt}{dx} (Q_i^{n} -  Q_{i-1}^{n})
		\end{equation*}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		\begin{figure}
			\includegraphics[width=0.3\textwidth]{\figpath/fig_fv_update.png}
		\end{figure}
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{Upwind scheme}
		\begin{block}{Stability criterion:}
		\begin{equation*}
			| \frac{a dt }{dx} | \leq \ 1
		\end{equation*}		
		\end{block}
		\vspace{1cm}
		Upwind scheme is of 1st order accuracy only and very dispersive. Therefore it is not accurate enough to be of any use for actual simulation tasks.
	\end{frame}	
				
% ---------------------------------------------------------

\subsection{The Lax-Wendroff Scheme}
	\begin{frame}{The Lax-Wendroff Scheme}
		Our goal is to find solutions to $\partial_t Q + a \partial_x Q = 0$. We start by using the Taylor expansion to extrapolate $Q(x,t)$ in time to obtain
		\begin{equation*}
			Q(x,t^{n+1}) = Q(x,t^n) + dt \partial_t Q(x,t^n) + \frac{1}{2} dt^2 \partial^2_t Q(x,t^n) + \dots 
		\end{equation*}
		From the governing equation we can also state by additional differentiations
		\begin{equation*}
			\begin{split}
				\partial^2_t Q &=  - a \partial_x \partial_t Q \\
				\partial_x \partial_t Q &=  \partial_t \partial_x Q = \partial_x (-a \partial_x Q) \\
				\partial^2_t Q &=  a^2 \partial^2_x Q
			\end{split}
		\end{equation*}
		noting that we just derived the 2nd order form of the acoustic wave equation.
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{The Lax-Wendroff Scheme}	
		Replacing time derivatives by the equivalent expressions containing space derivatives only and obtain
		\begin{equation*}
			Q(x,t^{n+1}) = Q(x,t^n) - dt a \partial_x Q(x,t^n) + \frac{1}{2} dt^2 a^2 \partial^2_x Q(x,t^n) + \dots
		\end{equation*}
		Using central differencing schemes for both space derivatives
		\begin{equation*}
			\begin{split}
				\partial_x Q(x,t^n) &\approx \frac{Q^n_{i+1} - Q^n_{i-1}}{2 dx} \\
				\partial^2_x Q(x,t^n) &\approx \frac{Q^n_{i+1} - 2 Q^n_i + Q^n_{i-1}}{dx^2}
			\end{split}
		\end{equation*}
		we finally obtain a fully discrete second-order scheme
		\begin{equation*}
			Q_i^{n+1} =   Q_i^{n} - \frac{a dt}{2 dx} (Q^n_{i+1} - Q^n_{i-1}) + \frac{1}{2} (\frac{a dt}{dx})^2 (Q^n_{i+1} - 2 Q^n_i + Q^n_{i-1})
		\end{equation*}
		known as the Lax-Wendroff scheme. 
	\end{frame}
				
% ---------------------------------------------------------

	\begin{frame}{The Lax-Wendroff Scheme}	
		\begin{itemize}
			\item[] Using standard finite-difference considerations without making use of flux concepts.
			\item[] By extending the finite-volume method towards higher order to approximate the solution inside the finite volume as a piecewise linear function.
			\item[] The choice of slope considered then determines the specific 2nd order numerical scheme that evolves.
		\end{itemize}
		The Lax-Wendroff scheme can also be interpreted as a finite-volume method by considering the flux functions
		\begin{equation*}
			\begin{split}
				F^n_L &= \frac{1}{2} a (Q^n_{i-1} + Q^n_{i}) -  \frac{1}{2}  \frac{dt}{dx} a^2 (Q^n_{i} - Q^n_{i-1})  \\ 
				F^n_R &= \frac{1}{2} a (Q^n_{i} + Q^n_{i+1}) -  \frac{1}{2}  \frac{dt}{dx} a^2 (Q^n_{i+1} - Q^n_{i}) 
			\end{split}
		\end{equation*}
	\end{frame}	
	
% ---------------------------------------------------------

\section{The Finite-Volume Method: Scalar Advection}

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

	\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		We proceed with implementing the two numerical schemes 1) the upwind method and 2) the Lax-Wendroff scheme. Recalling their formulations
		\begin{equation*}
			Q_i^{n+1} \ = \ Q_i^{n} - \frac{a dt}{dx} (Q_i^{n} -  Q_{i-1}^{n})	
		\end{equation*}
		and
		\begin{equation*}
			Q_i^{n+1} =   Q_i^{n} - \frac{a dt}{2 dx} (Q^n_{i+1} - Q^n_{i-1}) + \frac{1}{2} (\frac{a dt}{dx})^2 (Q^n_{i+1} - 2 Q^n_i + Q^n_{i-1})
		\end{equation*}
		Using a spatial initial condition, a Gauss function 
		\begin{equation*}
			Q(x,t=0) = e^{-1/\sigma^2 (x-x_0)^2}
		\end{equation*}
		that is advected with speed $c=2500 m/s$. The analytical solution to this problem is a simple translation of the initial condition to  $x = x_0+c t$, where $t=i dt$ is the simulation time at time step $i$.  	
	\end{frame}	
	
	
	% ---------------------------------------------------------

\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		\begin{minipage}{0.4\textwidth}
			\begin{table}
\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
$x_{max}$ & 75000 m \\
nx & 6000 \\
c & 2500 m/s \\
dt & 0.0025 s \\
dx & 12.5 m \\
$\epsilon$ & 0.9 \\
$\sigma$ (Gauss) & 200 m \\
$x_0$ & 1000 m
\end{tabular}
			\end{table}
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
	   \begin{figure}
			\includegraphics[width=0.9\textwidth]{./figs/fv_1.png}
		\end{figure}
		\end{minipage}		
	\end{frame}
	
% ------------------------------------

	\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		Implemement periodic and absorbing boundary conditions with the statements
		\begin{equation*}
			\begin{split}
				\text{Periodic: } \ Q^n_{1} \ &= Q^n_{nx-1} \\
				\text{Absorbing:} \ Q^n_{nx} \ &= Q^n_{nx-1} \  
			\end{split}
		\end{equation*}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_fv_abs.png}
			\caption{Boundary conditions. Absorbing, or circular boundary conditions can be implemented by using ghost cells outside the physical domain $x\in[x_0, x_{max}]$.}
		\end{figure}
	\end{frame}
					
% -------------------------------------------

	\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_fv_scalar.png}
			\caption{\footnotesize{{\bf Top:} Snapshots of an advected Gauss function (analytical solution, thick solid line) are compared with the numerical solution of the 1st order upwind method (thin solid line) and the 2nd order Lax-Wendroff scheme (dotted line) for increasing propagation distances. {\bf Bottom:} The same for a box-car function.}}
		\end{figure}
	\end{frame}	
					
% ---------------------------------------------------------

	\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		Implemement periodic and absorbing boundary conditions with the statements
		\begin{equation*}
			\begin{split}
				\text{Periodic: } \ Q^n_{1} \ &= Q^n_{nx-1} \\
				\text{Absorbing:} \ Q^n_{nx} \ &= Q^n_{nx-1} \  
			\end{split}
		\end{equation*}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_fv_abs.png}
			\caption{Boundary conditions. Absorbing, or circular boundary conditions can be implemented by using ghost cells outside the physical domain $x\in[x_0, x_{max}]$.}
		\end{figure}
	\end{frame}
					
% ---------------------------------------------------------

	\begin{frame}{The Finite-Volume Method: Scalar Advection}	
		\begin{figure}
			\includegraphics[width=0.6\textwidth]{\figpath/fig_fv_scalar.png}
			\caption{\footnotesize{{\bf Top:} Snapshots of an advected Gauss function (analytical solution, thick solid line) are compared with the numerical solution of the 1st order upwind method (thin solid line) and the 2nd order Lax-Wendroff scheme (dotted line) for increasing propagation distances. {\bf Bottom:} The same for a box-car function.}}
		\end{figure}
	\end{frame}	
					
% ---------------------------------------------------------	

\section{The Finite-Volume Method for Elastic Waves}

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

	\begin{frame}{The Finite-Volume Method for Elastic Waves}
		\begin{block}{Source-free version of the coupled first-order elastic wave equation}
			\begin{equation*}
				\begin{split}
					\partial_t v - \frac{1}{\rho} \partial_x \sigma \ &= \ 0 \\  
					\partial_t {\sigma} \ - \ \mu \partial_x v \ &= \ 0  \ .
				\end{split}
			\end{equation*}
		\end{block}
		\vspace{0.5cm}
		We proceed by writing this equation in matrix-vector notation
		\begin{equation*}
			\partial_t \mathbf{Q} + \mathbf{A} \partial_x \mathbf{Q} = 0 
		\end{equation*}
		where $\mathbf{Q}=(\sigma, v)$ is the vector of unknowns and matrix $\mathbf{A}$ contains the parameters
		\begin{equation*}
			\mathbf{A} \ = \
			\begin{pmatrix}
				0 & -1/\rho \\
				-\mu & 0 
			\end{pmatrix}
		\end{equation*}
	\end{frame}		
						
% ---------------------------------------------------------	

	\begin{frame}{The Finite-Volume Method for Elastic Waves}
		The problem hereby is that these eqations are coupled. What needs to be done is to demonstrate the {\it hyperbolicity} of the wave equation in this form, i.e. show that $\mathbf{A}$ is diagonalizable.\\
		In the case of a quadratic matrix $\mathbf{A}$ with shape $m \times m$ 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 Finite-Volume Method for Elastic Waves}
		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 above equation we obtain 
		\begin{equation*}
			\textbf{R}^{-1}\partial_t \mathbf{Q} + \textbf{R}^{-1}\textbf{R}\Lambda\textbf{R}^{-1} \partial_x \mathbf{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 Finite-Volume Method for Elastic Waves}
		What remains to be shown is that in our specific case \textbf{A} has real eigenvalues. These are easily determined as $\lambda_{1,2}=\pm \sqrt{\mu/\rho}= \pm c$, corresponding to the shear velocity $c$. For the eigenvectors we obtain
		\begin{equation*}
			\textbf{r}_{1,2} \ = \
			\begin{pmatrix}
			\pm \rho c \\
			1
			\end{pmatrix}
		\end{equation*}
		interestingly enough containing as first elements values of the {\it seismic impedance} $Z=\rho c$ that is 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 Finite-Volume Method for Elastic Waves}		
		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.\\
		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 Finite-Volume Method for Elastic Waves}		
		Obtaining the final analytical solution for velocity $v$ and stress $\sigma$ 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 this solution can be expressed as  
		\begin{equation*} 
			\textbf{Q}(x,t) \ = \ \sum_{p=1}^m \mathbf{w}_p(x,t) \mathbf{r}_p 
		\end{equation*}
		meaning that any solution is a sum over weighted eigenvectors, a superposition of $m$ waves. 
	\end{frame}
						
% ---------------------------------------------------------	

\subsection{Homogeneous Case}
	\begin{frame}{Homogeneous Case}
		\begin{minipage}{0.5\textwidth}
			\begin{figure}
				\includegraphics[width=0.9\textwidth]{\figpath/fig_fv_riemann_homo.png}		
			\end{figure}
		\end{minipage}		
    \mbox{ }
    \hfill	
		\begin{minipage}{0.45\textwidth}
			\begin{block}{}
				\footnotesize{Riemann problem, homogeneous case. {\bf Top:} A discontinuity  $\Delta Q$ is located at $x=0$ as initial condition to the advection equation (e.g., as initial stress discontinuity). {\bf Bottom:} The discontinuity propagates along characteristic curves in the space-time domain. The figure illustrates adjacent cells and two time levels $t_n$ and $t_{n+1}$. Two waves propagate in opposite direction modifying the values in the cells adjacent to $x=0$.}
			\end{block}
		\end{minipage}	
	\end{frame}	
						
% ---------------------------------------------------------	

	\begin{frame}{Homogeneous Case}
		The solution to our problem is a superposition of weighted eigenvectors $\mathbf{r}_p$, in our case $p=1, 2$. Therefore, we can decompose the discontinuity jump into these eigenvectors according to
	\begin{equation*}
		\begin{split}
			\Delta \mathbf{Q} = \mathbf{Q}_r -\mathbf{Q}_l &= \alpha_1 \mathbf{r}_1 + \alpha_2 \mathbf{r}_2  \\
			  \mathbf{R} \alpha & =  \Delta \mathbf{Q} \\
			  \alpha & =  \mathbf{R}^{-1} \Delta \mathbf{Q} 
		\end{split}
	\end{equation*}
	where $\mathbf{R}$ is the matrix of eigenvectors as defined above and $\alpha$ are unknown weights.
	\end{frame}
						
% ---------------------------------------------------------	

	\begin{frame}{Homogeneous Case}
		Decompose the solution into positive (right-propagating) and negative (left-propagating) eigenvalues
		\begin{equation*}
			\Lambda^- =
			\begin{pmatrix}
			-c & 0 \\
			0 & 0	
			\end{pmatrix}
			\ , \ \ 
			\Lambda^+ =
			\begin{pmatrix}
			0 & 0 \\
			0 & c	
			\end{pmatrix}  
		\end{equation*}
		Then we can derive matrices $\mathbf{A}^\pm$ - corresponding to the advection velocity in the scalar case  
		\begin{equation*}
			\begin{split}
				\mathbf{A}^+ &= \mathbf{R} \Lambda^+ \mathbf{R}^{-1} \\
				\mathbf{A}^- &= \mathbf{R} \Lambda^- \mathbf{R}^{-1}
			\end{split}
		\end{equation*}
	\end{frame}	
						
% ---------------------------------------------------------	

	\begin{frame}{Homogeneous Case}		
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_fv_flux_homo.png}
			\caption{\footnotesize{The constant average cell values of $Q$ are illustrated for three adjacent cells $i-1, i, i+1$. The eigenvector decomposition leads to wave $A^+$ propagating from the left boundary with velocity $c$ in to cell $i$ and wave $A^-$ propgating with velocity $-c$ into cell $i$ from the right boundary. This determines the flux of discontinuities $\Delta \mathbf{Q}_{l,r}$ into cell $i$ by the amount $dt/dx \Delta \mathbf{Q}_{l,r}$.}}
		\end{figure}
	\end{frame}		
						
% ---------------------------------------------------------	

	\begin{frame}{Homogeneous Case}	
		We are  ready to formulate an upwind finite-volume scheme for any multi-dimensional linear hyperbolic system:
		\begin{equation*}
			\begin{split}
				\Delta \mathbf{Q}_l &= \mathbf{Q}_i - \mathbf{Q}_{i-1} \\
				\Delta \mathbf{Q}_r &= \mathbf{Q}_{i+1} - \mathbf{Q}_i \\
				\mathbf{Q}_i^{n+1} &= \mathbf{Q}^n_i - \frac{dt}{dx} (\mathbf{A}^+ \Delta \mathbf{Q}_l + \mathbf{A}^- \Delta \mathbf{Q}_r) \ . 
			\end{split}
		\end{equation*}
		We can relate this formulation to the very basic flux concept
		\begin{equation*}
			\begin{split}
				\mathbf{F}_{l} \ &= \ \mathbf{A}^+ \Delta \mathbf{Q}_l\\ 
				\mathbf{F}_{r} \ &= \ \mathbf{A}^- \Delta \mathbf{Q}_r \ . 
			\end{split}
		\end{equation*}
	\end{frame}
							
% ---------------------------------------------------------	

	\begin{frame}{Homogeneous Case}	
		The 1st order upwind solution is of no practical use because of its strong dispersive behaviour. \\
		$\Longrightarrow$The 2nd order Lax-Wendroff scheme\\
		\vspace{0.5cm}
		The high-order scheme does not necessitate the separation into eigenvectors and the Jacobian matrix $A$ can be used in its original form. The extrapolation scheme reads
		\begin{equation*}
			\begin{split}
				\mathbf{Q}_i^{n+1} =& \mathbf{Q}^n_i - \frac{dt}{2 dx} \mathbf{A} (\mathbf{Q}^n_{i+1} +  \mathbf{Q}^n_{i-1}) \\
				&\ \ \ \ + \frac{1}{2}\frac{dt}{dx}^2 \mathbf{A}^2 (\mathbf{Q}^n_{i-1} - 2 \mathbf{Q}_i^n +  \mathbf{Q}^n_{i+1}) \ . 
			\end{split}
		\end{equation*}
	\end{frame}
	
	
% ---------------------------------------------------------	

	\begin{frame}{Example}
				\begin{minipage}{0.4\textwidth}
			\begin{table}
\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
$x_{max}$ & 10000 m \\
nx & 800 \\
c & 2500 m/s \\
$\rho$  & 2500 $kg/m^3$ \\
dt & 0.025 s \\
dx & 12.5 m \\
$\epsilon$ & 0.5 \\
$\sigma$ (Gauss) & 200 m \\
$x_0$ & 5000 m
\end{tabular}
			\end{table}
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=1.\textwidth]{./figs/fv_2.png}
		\end{figure}	
		\end{minipage}		
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Result}
		\begin{minipage}{0.4\textwidth}		
		\begin{figure}
			\includegraphics[width=0.55\textwidth]{\figpath/fig_fv_ls_homo.png}
		\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
			\begin{block}{}
				The stress-velocity system is advected for an initial condition of Gaussian shape (top, dashed line, scaled by factor 1/2). \textbf{Top: }Stress snapshot at time $t=1.5 s$.  \textbf{Bottom:} Velocity snapshot at the same time. In both cases analytical solutions are superimposed.
			\end{block}
		\end{minipage}
	\end{frame}	
							
% ---------------------------------------------------------	

\subsection{Heterogeneous Case}
	\begin{frame}{Heterogeneous Case}
		\begin{minipage}{0.4\textwidth}		
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fv_riemann_hetero.png}
		\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
			\begin{block}{}
				{\bf Top:} Single discontinuity separating two regions with different properties. {\bf Bottom:} Velocities and impedances on both sides of the discontinuity. The Riemann problem solves the problem of how waves on both sides are partitioned.
			\end{block}
		\end{minipage}	
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}
		Mathematically the solution still has to consist of a weighted sum over eigenvectors that now describe solutions in the left and right parts.
		\begin{equation*}
			\Delta \mathbf{Q} = \alpha_1 
			\begin{pmatrix}
			-Z_l \\
			1
			\end{pmatrix}
			+
			\alpha_2 
			\begin{pmatrix}
			Z_r \\
			1
			\end{pmatrix}		
		\end{equation*}
		for some unknown scalar values $\alpha_{1,2}$. This can be written as a linear system of the  form
		\begin{equation*}
			\mathbf{R}_{lr} \alpha = \Delta \mathbf{Q}
		\end{equation*}
		where $\alpha$ is a vector and the matrices with the eigenvector are
		\begin{equation*}
			\mathbf{R}_{lr} = 
			\begin{pmatrix}
			-Z_l & Z_r \\
			1 & 1
			\end{pmatrix}
			\ , \ \ 
			\mathbf{R}^{-1}_{lr} = \frac{1}{Z_l + Z_r} 
			\begin{pmatrix}
			-Z_l & Z_r \\
			1 & 1
			\end{pmatrix}
		\end{equation*}
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}
		Using the $\mathbf{R}$ matrix for cell boundaries separating cells with different (constant) properties. Again we separate into left- and right-propagating eigenvalues within cell $i$ 
		\begin{equation*}
			\Lambda^- =
			\begin{pmatrix}
			-c_i & 0 \\
			0 & 0	
			\end{pmatrix}
			\ , \ \ 
			\Lambda^+ =
			\begin{pmatrix}
			0 & 0 \\
			0 & c_i	
			\end{pmatrix}  		
		\end{equation*}
		and using the definitions 
		\begin{equation*}
			\mathbf{R}_{l} = 
			\begin{pmatrix}
			-Z_{i-1} & Z_i \\
			1 & 1
			\end{pmatrix}
			\ , \ \ 
			\mathbf{R}_{r} = 
			\begin{pmatrix}
			-Z_{i} & Z_i+1 \\
			1 & 1
			\end{pmatrix}		
		\end{equation*}
		for the eigenvectors describing the solutions around the left and right boundaries we can determine the corresponding advection terms as
		\begin{equation*}
			\begin{split}
				\mathbf{A}^+ &= \mathbf{R}_l \Lambda_l^+ \mathbf{R}_l^{-1} \\
				\mathbf{A}^- &= \mathbf{R}_r \Lambda_r^+ \mathbf{R}_r^{-1}
			\end{split}	
		\end{equation*}
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}		
		Leading to the 1st order upwind  extrapolation scheme for the solution vector $Q_i$ in the general heterogeneous case
		\begin{equation*}
			\begin{split}
				\Delta \mathbf{Q}_l &= \mathbf{Q}_i - \mathbf{Q}_{i-1} \\
				\Delta \mathbf{Q}_r &= \mathbf{Q}_{i+1} - \mathbf{Q}_i \\
				\mathbf{Q}_i^{n+1} &= \mathbf{Q}^n_i - \frac{dt}{dx} (\mathbf{A}^+ \Delta \mathbf{Q}_l + \mathbf{A}^- \Delta \mathbf{Q}_r) 
			\end{split}
		\end{equation*}
		 $\Longrightarrow$ too dispersive but interesting side effect!\\
		 \vspace{0.3cm}
		  Let us take the eigenvector (i.e., wave) propagating in the left domain. What does this imply for the wave propagating in the right domain?
		\begin{equation*}
			\Delta \mathbf{Q} = 
			\begin{pmatrix}
			Z_l \\
			1
			\end{pmatrix}		
		\end{equation*}
	\end{frame}	
								
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}	
		We need to determine how the waves are partitioned given the discontinuous material parameters described by matrix $\mathbf{R}_{lr}$. This leads to the coefficients $\alpha$ as
		\begin{equation*}
			\begin{split}
				\alpha &= \mathbf{R}^{-1}_{lr} \Delta \mathbf{Q} \\
				&= \frac{1}{Z_l+Z_r}
				\begin{pmatrix}
				-1 & Z_r \\
				1 & Z_l
				\end{pmatrix}
				\begin{pmatrix}
				Z_l \\
				1 
				\end{pmatrix} 
				\\
				&=
				\renewcommand\arraystretch{2}
				\begin{pmatrix} 
				\frac{Z_r-Z_l}{Z_l+Z_r}\\
				\frac{2 Z_l}{Z_l+Z_r}
				\end{pmatrix}   
				=
				\begin{pmatrix} 
				R \\
				T
				\end{pmatrix} 
			\end{split}		
		\end{equation*}
		that are the well known transmission and reflection coefficients for vertical incidence at a material discontinuity.
	\end{frame}
									
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}	
		\begin{minipage}{0.4\textwidth}		
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fv_refcoef.png}
		\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
			\begin{block}{}
				\textbf{Reflection and transmission coefficients}. Seismic waves incident perpendicular to a material discontinuity are reflected and transmitted according to coefficients $R$ and $T$. These coefficients can be derived via the Riemann problem used to develop flux schemes for finite-volume methods.
			\end{block}
		\end{minipage}		
	\end{frame}
									
% ---------------------------------------------------------	

	\begin{frame}{Heterogeneous Case}
		Finally, we present the 2nd order scheme for the general heterogeneous results and a simulation example.
		\begin{equation*}
			\begin{split}
				\mathbf{A}_{i-1/2} &= \frac{1}{Z_{i-1}+Z_i} 
				\begin{pmatrix}
				c_i Z_i - c_{i-1} Z_{i-1}  & (c_{i-1} + c_i) Z_{i-1} Z_{i} \\
				c_ {i-1}  + c_{i}  & c_{i} Z_{i-1} - c_{i-1} Z_{i} 
				\end{pmatrix}
				\\
				\mathbf{A}_{i+1/2} &= \frac{1}{Z_{i}+Z_{i+1}} 
				\begin{pmatrix}
				c_{i+1} Z_{i+1} - c_{i} Z_{i}  & (c_{i} + c_{i+1}) Z_{i} Z_{i+1} \\
				c_i  + c_{i+1}  & c_{i+1} Z_i - c_i Z_{i+1} 
				\end{pmatrix}  \ .
			\end{split}
		\end{equation*}
		With these definitions we can formulate the Lax-Wendroff extrapolation scheme for elastic waves in heterogeneous material 
		\begin{equation*}
			\begin{split}
				\mathbf{Q}_i^{n+1} =& \mathbf{Q}^n_i - \frac{dt}{2 dx} \left[ \mathbf{A}_{i-1/2} (\mathbf{Q}^n_{i} -  \mathbf{Q}^n_{i-1} ) +   											\mathbf{A}_{i+1/2} (\mathbf{Q}^n_{i+1} -  \mathbf{Q}^n_{i} ) \right] \\
				& \ \ \ \ + \frac{1}{2}(\frac{dt}{dx})^2 \left[ \mathbf{A}_{i+1/2}^2 (\mathbf{Q}^n_{i+1} - \mathbf{Q}^n_i) -  \mathbf{A}_{i-1/2}^2 											(\mathbf{Q}^n_{i} - \mathbf{Q}^n_{i-1})\right] 
			\end{split}
		\end{equation*}
	\end{frame}	
										
% ---------------------------------------------------------	

	\begin{frame}{Example}
				\begin{minipage}{0.38\textwidth}
			\begin{table}
\begin{tabular}[t]{ll}
Parameter & Value \\
\hline
$x_{max}$ & 10000 m \\
nx & 800 \\
$c_l$ & 2500 m/s \\
$c_r$ & 5000 m/s\\
$\rho$  & 2500 $kg/m^3$ \\
dt & 0.025 s \\
dx & 12.5 m \\
$\epsilon$ & 0.5 \\
$\sigma$ (Gauss) & 200 m \\
$x_0$ & 5000 m
\end{tabular}
			\end{table}
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.57\textwidth}
		\begin{figure}
			\includegraphics[width=1\textwidth]{./figs/fv_3.png}
		\end{figure}	
		\end{minipage}				
	\end{frame}
										
% ---------------------------------------------------------	

	\begin{frame}{Result}
		\begin{figure}
			\includegraphics[width=0.55\textwidth]{\figpath/fig_fv_ls_hetero.png}
		\end{figure}
	\end{frame}	
										
% ---------------------------------------------------------	
	
\section{The Finite-Volume Method via Gauss' Theorem}	
											
% ---------------------------------------------------------	

\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
Gauss' theorem states that the outward flux of a vector field $Q_i(x,t)$ through a closed surface $S$ is equal to the volume integral of the divergence over the volume $V$ inside the surface at some time $t$. Mathematically this can be expressed as 
\begin{equation*}
\int_V \partial_i Q_i dV = \int_S n_i Q_i dS
\end{equation*}
where $n_i$ are the components of the local surface normal vector. \\
Note:
\begin{itemize}
	\item[] Not restricted to vector fields $\longrightarrow$ any tensor field $Q_{i,j\ldots}$
	\item[] Applies also to scalar fields Q(x,t)
\end{itemize}
\end{frame}
										
% ---------------------------------------------------------	

\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
Assuming the gradient of the solution field is smooth enough and can be assumed constant inside volume $V$ we can take it out of the integral and obtain an expression for the derivative as a function of an integral over a surface $S$ with segements $dS$ in 3D or a line with segments $dL$ surrounding a surface $S$ in 2D
\begin{equation*}
\begin{split}
\partial_i Q \int_V dV &= \int_S n_i Q dS \\
\partial_i Q_{3D} &= \frac{1}{V} \int_S n_i Q dS \\
\partial_i Q_{2D} &= \frac{1}{S} \int_L n_i Q dL
\end{split}
\end{equation*}
\end{frame}

% ---------------------------------------------------------	
	
	\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
		\begin{minipage}{0.4\textwidth}		
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fv_gauss.png}
		\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
				Concept of numerical gradient calculation using Gauss' theorem illustrated in 2D. The constant gradient of a scalar field Q inside the {\it finite volume} $S$ is approximated as $\partial_i Q=1/S \sum_{\alpha} n_i^{\alpha} dL^{\alpha} Q^{\alpha}$. The polygon can have any shape.
		\end{minipage}			
	\end{frame}		
		
% ---------------------------------------------------------	

	\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
		\begin{minipage}{0.4\textwidth}		
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fv_lozenge.png}
		\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
Once a discretization of surfaces or surrounding lines is found we can develop a discrete scheme replacing the integrals in the above equations by sums to obtain  
\begin{equation*}
\begin{split}
\partial_i Q_{3D} &\approx \frac{1}{S} \sum_{\alpha} n_i^{\alpha} dS^{\alpha} Q^{\alpha} \\
\partial_i Q_{2D} &\approx \frac{1}{L} \sum_{\alpha} n_i^{\alpha} dL^{\alpha} Q^{\alpha}  \  . 
\end{split}
\end{equation*}
The importance here is that this description is entirely independent of the shape of a particular volume.
		\end{minipage}			
	\end{frame}	
		
% ---------------------------------------------------------	

\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
A rhombus-shaped finite volume cell is described by four edge points.  
\begin{equation*}
\begin{split}
P_1 = (-\Delta_1,0), P_2 = (+\Delta_1,0), \\
P_3 = (0,-\Delta_2), P_4 = (0,+\Delta_2)
\end{split}
\end{equation*}
with the length of the sides given by $\ell=\sqrt{\Delta^2_1+\Delta^2_2}$ and the surface $S=2\Delta_1\Delta_2$. The four normal vectors are defined by
\begin{equation*}
\begin{split}
n_{2,4} = 
\begin{pmatrix}
\Delta_2 \\
\Delta_1
\end{pmatrix}; \ 
&n_{4,1} = 
\begin{pmatrix}
-\Delta_2 \\
\Delta_1
\end{pmatrix}; \\
n_{1,3} = 
\begin{pmatrix}
-\Delta_2 \\
-\Delta_1
\end{pmatrix}; \ 
&n_{3,2} = 
\begin{pmatrix}
\Delta_2 \\
-\Delta_1
\end{pmatrix} \ .
\end{split}
\end{equation*}
\end{frame}
		
% ---------------------------------------------------------	

\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
We now have all components to apply previous equation in the 2D case. We integrate along the paths suggested in the figure to obtain 
\begin{equation*}
\begin{split} 
\partial_1 Q &= \frac{1}{S}(  \frac{\ell}{2} Q_1 ( -  \frac{\Delta_2}{\ell} - \frac{\Delta_2}{\ell}) + \frac{\ell}{2} Q_2(\frac{\Delta_2}{\ell}+ \frac{\Delta_2}{\ell}) \\
&\ + \frac{\ell}{2} Q_3 ( -  \frac{\Delta_2}{\ell} + \frac{\Delta_2}{\ell}) + \frac{\ell}{2} Q_4(\frac{\Delta_2}{\ell}- \frac{\Delta_2}{\ell})) \\
&=\frac{1}{S}(-Q_1\Delta_2+Q_2\Delta_2) \\
&=\frac{Q_2-Q_1}{2\Delta_1} \ 
\end{split}
\end{equation*}
for the 1st derivative of $Q$ w.r.t. $x$. 
\end{frame}
		
% ---------------------------------------------------------	

\begin{frame}{The Finite-Volume Method via Gauss' Theorem}
\begin{minipage}{0.4\textwidth}
	\begin{figure}
		\includegraphics[width=0.65\textwidth]{\figpath/fig_fv_kaeser.png}
	\end{figure}
		\end{minipage}
    \mbox{ }
    \hfill			
		\begin{minipage}{0.55\textwidth}
Example of elastic wave simulations in 2D using difference operators based on the finite-volume approach (Kaeser 2001b). In this example an unstructured grid follows a free-surface topography with ghost cells outside the surface to implement stress-free boundary conditions.
		\end{minipage}
\end{frame}

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

\begin{frame}{Summary}
	\footnotesize{
	\begin{itemize}
		\item The finite-volume method naturally follows from discretizing conservation equations considering fluxes between finite-volume cells of averaged solution fields.
\item The fluxes across boundaries during an extrapolation step are estimated using solutions to the Riemann problem.
\item The Riemann problem considers the advection of a single jump discontinuity taking into account the analytical solution of the homogeneous problem. It allows an analytical prediction of how much of the material (energy, stress, etc.) to be conserved enters into or leaves a cell.
\item The lowest order finite-volume solution to the advection equation leads to a finite-difference algorithm with a forward (or backward) spatial differencing scheme, depending on the advection direction. This is called an upwind scheme.
	\end{itemize}
	}
\end{frame}

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

\begin{frame}{Summary}
	\footnotesize{
	\begin{itemize}
\item First-order finite-volume schemes are highly dispersive and not appropriate for the solution of wave propagation problems. The 2nd order Lax-Wendroff scheme does a much better job.
\item The problem of elastic wave propagation can be formally cast as a 1st order hyperbolic problem. Therefore - only with slight modifications - the fundamental schemes developed for the scalar advection problem can be applied to elastic wave propagation.
\item In the finite-volume method the problem of estimating partial derivatives (finite differences) is replaced by the requirement to accurately calculate fluxes across cell boundaries.
\item A main advantage of the finite-volume method is the fact that the scheme can be easily applied to volume cells of any shape.
\item Finite-volume schemes for arbitrary high-order reconstructions inside the cells and high-order time-extrapolation schemes have been developed but not used extensively in seismology.
	\end{itemize}
	}
\end{frame}	
		
% ---------------------------------------------------------	

\begin{frame}{Comprehension questions}
\footnotesize{
	\begin{enumerate}
\item[1] What is the connection between finite-volume methods and conservation equations?
\item[2] What is meant by finite {\it volumes}, is there any difference to a finite {\it element}?
\item[3] If you look at the uwpind approach to the scalar advection problem (Eq.~\ref{fv_discrete}), why is the finite-volume method so closely linked to staggered-grid finite-difference schemes? Explain.
\item[4] What are the main  advantages of finite-volume methods compared with finite-difference methods?
\item[5] Explain in a qualitative way what the Riemann problem is and why it is so essential for finite-volume schemes. 
\item[6] In what areas of natural sciences are finite-volume schemes mostly used. Explore the literature and try to give reasons.
\item[7] What is numerical diffusion? Why is it relevant for finite-volume methods?
\item[8] What is the connection between reflection/transmission coefficients of seismic waves and the finite-volume method?
\item[9] The finite-volume method extrapolates cell averages. What strategies do you see to extend the method to high-order accuracy?
\end{enumerate}
}
\end{frame}	
		
% ---------------------------------------------------------	

\begin{frame}{Theoretical problems}
\footnotesize{
\begin{enumerate}
\item[10] Show that equation is a finite-difference solution to the equation $\partial_t Q-a\partial_x Q=0$ using a forward difference in space. 
\item[11] The stability criterion for the finite-volume method is $c dt/dx \leq 1$.  Starting with Fig. derive this stability criterion from first principles. 
\item[12] Starting with the advection equation $\partial_t Q-a\partial_x Q=0$ derive the second-order wave equation by applying the so-called Cauchy-Kovalevskaya procedure.
\item[13] Following the finite-volume approach based on the divergence theorem calculate the spatial derivative operator for the hexagonal cell.
\end{enumerate}
}
\end{frame}	
		
% ---------------------------------------------------------	

\begin{frame}{Theoretical problems}
\footnotesize{
\begin{enumerate}
\item[14] The linear system for elastic wave propagation in 1D (transverse motion) is given in equation. The wave equation can also be formulated for compressional waves using the compressibility $K$ as elastic constant. Reformulate the linear system for acoustic wave propagation and calculate the eigenvalues of the resulting Jacobian matrix $\mathbf{A}$.
\item[15] For either elastic or acoustic linear system derive the eigenvectors of Jacobian matrix $\mathbf{A}$, the matrix of eigenvectors and its inverse. 
\item[16] Show that the superposition of left- and right-propagating stress and velocity waves are solutions to the linear system of equations for elastic wave propagation. 
\item[17] Show that a discontinuity of the form $\Delta Q = [1,0]$ leads to an equi-partitioning of two seismic waves propagating in opposite directions. Start with the Riemann problem formulated for the homogeneous case.
\item[18] Derive reflection and transmission coefficients for seismic waves with vertical incidence by considering the Riemann problem for material discontinuity.
\end{enumerate}
}
\end{frame}	
			
% ---------------------------------------------------------	

\begin{frame}{Programming Exercises}
\footnotesize{
\begin{enumerate}
\item[19] Determine the stability limit of the Lax-Wendroff scheme for the scalar advection equation. Is it the same as for the upwind scheme?
\item[20] Create a highly unstructured 1D mesh and investigate the accuracy of the finite-volume method (Lax-Wendroff) for the scalar advection problem.
\item[21] Investigate the concept of trapped elastic waves by injecting an initial condition in a low-velocity region. Use the Lax-Wendroff algorithm in 1D.  
\item[22] Implement circular boundary conditions in the 1D elastic Lax Wendroff solution. Initiate a sinusoidal function $f(x)=sin(kx)$ that is advected in one direction. Investigate the accuracy of the finite-volume scheme as a function of wavelength and propagation distance by comparing with the analytical solution.   
\item[23] The finite-volume method is supposed to conserve energy in the homogeneous case. Use the computer programs for scalar advection, set up an example and calculate the total energy in the system for each time step. Check whether it is conserved. Explore this problem for the heterogeneous case.   
\item[24] Scalar advection problem: Advect a Gaussian shaped waveform as long as you can and extract the travel time difference with the analytical solution in an automated way using cross-correlation. Plot the time error as a function or propagation distance and your simulation parameters.
\end{enumerate}
}
\end{frame}	


\end{document}
