% beamer version of chapter 4 "The Finite Difference Method" in "Computational Seismology, A Practical Introduction" (H. Igel)
%
%\documentclass{beamer}
%
%\mode<presentation>
%{
  %\usetheme{CambridgeUS}
%
  %\setbeamercovered{transparent}
  %\setbeamertemplate{navigation symbols}{}
%}

\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{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_fd}
\newcommand{\fs}{\footnotesize}

\title[Computational Seismology] % (optional, use only with long paper titles)
{The Finite Difference 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{Stability, Dispersion, Convergence}
\subsection{Stability}

	\begin{frame}{Stability}
		Starting with the definition of a plane harmonic wave for pressure p propagating in x-direction with wavenumber k and angular frequency $\omega$
		\begin{equation*}
			p(x, t) \ = \ e^{i(kx - \omega t)}
		\end{equation*}
		In the finite-difference approximation to the wave equation, the spatio-temporal discretization was
		\begin{equation*}
			\begin{split}
				x \ &\rightarrow \ j \mathrm{d}x \\
				t \ &\rightarrow \ n \mathrm{d}t 
			\end{split}
		\end{equation*}
	\end{frame}

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

	\begin{frame}{Stability}
		Using this for the plane-wave formula such that
		\begin{equation*}
			\begin{split}
				P_j^n \ &\rightarrow \ e^{i(k j \mathrm{d} x - \omega n \mathrm{d}t)} \\
				P_{j+1}^n \ &= \ e^{i(k (j+1) \mathrm{d} x - \omega n \mathrm{d}t)} \\
        							  \ &= \ e^{i k \mathrm{d} x} e^{i(k j \mathrm{d} x - \omega n \mathrm{d}t)} \\
        							  \ &= \ e^{i k \mathrm{d} x} P_j^n \\
				P_j^{n+1} \ &= \ e^{-i \omega \mathrm{d} t} P_j^n
			\end{split}
		\end{equation*}
	\end{frame}

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

	\begin{frame}{Stability}
		Obtaining from the (source-free) finite-difference approximation of the acoustic wave equation and discretized plane-wave Ansatz
		\begin{equation*}
		e^{i \omega \mathrm{d}t} + e^{-i \omega  \mathrm{d}t} -2 \ = \ c^2 \frac{\mathrm{d}t^2}{\mathrm{d}x^2}
    	(e^{i k \mathrm{d}x} + e^{-i k \mathrm{d}x} -2 )
		\end{equation*}
		Dividing both sides by the term $P_j^n$ and using the definition
		\begin{equation*}
			\cos{x} \ = \ \frac{1}{2} (e^{i x} + e^{-i x} )
		\end{equation*}
		yields to
		\begin{equation*}
			\cos{( \omega \mathrm{d} t)} - 1 \ = \ c^2 \frac{\mathrm{d}t^2}{\mathrm{d}x^2} (\cos(k \mathrm{d} x) - 1)
		\end{equation*}
	\end{frame}

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

	\begin{frame}{Stability}
		and with the trigonometric relation
		\begin{equation*}
			\sin{ \frac{x}{2} } \ = \ \pm \sqrt{\frac{1 - \cos{x}}{2}}
		\end{equation*}
		finally arrive at
		\begin{block}{}
			\begin{equation*}
				\sin{( \omega \frac{\mathrm{d}t}{2} )} \ = \ c \frac{\mathrm{d}t}{\mathrm{d}x} \sin{(k \frac{\mathrm{d}x}{2} ) }
			\end{equation*}		
		\end{block}
		The equation only has real solutions when the connecting term has the following property
			\begin{equation*}
				\qquad  \boxed{c \frac{\mathrm{d}t}{\mathrm{d}x} \ \leq  \ 1 }
			\end{equation*}
	\end{frame}

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

	\begin{frame}{Courant-Friederichs-Lewy (CFL-) criterion}
		\begin{itemize}[<+->]
			\item The space-time discretization can not be arbitrarily chosen but depends on the medium properties (here: velocity $c$).
			\item The CFL-criterion describes the {\it conditional stability} that  allows convergent behaviour of the solution. 
			\item As the space discetization is often imposed by considering the smallest seismic velocities in the medium and the highest frequencies, the CFL criterion determines the time increment and thus the number of time steps to achieve a certain simulation length. 
			\item The actual value that has to be respected depends on the number of space dimensions and the overall algorithm.
			\item It is important to note that the fulfillment of the CFL criterion is by no means a guarantee for an accurate simulation!
		\end{itemize}
	\end{frame}

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

\subsection{Numerical Dispersion}

	\begin{frame}{Numerical Dispersion}
		The angular frequency $\omega$ can be expressed as
		\begin{equation*}
			\omega \ = \ \frac{2}{\mathrm{d}t} \sin^{-1}{ \left[ c \frac{\mathrm{d}t}{\mathrm{d}x} \sin{k \frac{\mathrm{d}x}{2}} \right] }
		\end{equation*}
		In this case the phase velocity should be identical to the acoustic velocity obtained by relating it to the wavenumber k
		\begin{equation*}
			c(k) \ = \ \frac{\omega}{k} \ = \ \frac{2}{k \mathrm{d}t} \sin^{-1}{ \left[ c \frac{\mathrm{d}t}{\mathrm{d}x} 
    \sin{k \frac{\mathrm{d}x}{2}} \right] }
		\end{equation*}
	\end{frame}

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

	\begin{frame}{Numerical Dispersion}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_courant_eps.png}
		\end{figure}
	\end{frame}

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

\subsection{Convergence}

	\begin{frame}{Convergence}
		When the finite-difference becomes infitesimally small the analytical derivative is recovered. Does that hold for our numerical approximation of the acoustic wave equation?
		\begin{equation*}
			c(k) \ = \ \frac{\omega}{k} \ = \ \frac{2}{k \mathrm{d}t} \sin^{-1}{ \left[ c \frac{\mathrm{d}t}{\mathrm{d}x} 
    \sin{k \frac{\mathrm{d}x}{2}} \right] }
		\end{equation*}
		\begin{equation*}
			\begin{split}
				\sin{x} \ &\approx \ x \qquad \text{for small} \ x \\
				\sin^{-1}{x} \ &\approx \ x \qquad \text{for small} \ x   \ . 
			\end{split}
		\end{equation*}
		As the spatial increment $dx$ and the time increment $dt$ converge to zero the original dispersion relation is recovered
		\begin{equation*}
			\lim_{\substack{\mathrm{d}t \to 0\\\mathrm{d}x \to 0}} c( \omega) \ = \ \frac{\omega}{k} \ = \ c_{\text{exact}}
		\end{equation*}
	\end{frame}
	
% ---------------------------------------------------------

\section{Acoustic Wave Propagation in 2D}


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

	\begin{frame}{Acoustic Wave Propagation in 2D}
		In 2D the constant-density acoustic wave equation is given as
		\begin{block}{}
			\begin{equation*}
				\ddot{p}(x,z,t) \ = \ c(x,z)^2 (\partial_x^2 p(x,z,t) + \partial_z^2 p(x,z,t)) \ + s(x,z,t)
			\end{equation*}
		\end{block}
		\begin{figure}
			\includegraphics[width=0.3\textwidth]{\figpath/fig_grid2d.png}
		\end{figure}
	\end{frame}	
	
% ---------------------------------------------------------	

	\begin{frame}{Acoustic Wave Propagation in 2D}	
		Discretizing space-time using
		\begin{equation*}
			p(x,z,t) \ \rightarrow \ p^n_{j,k} \ = \ p(n dt, j dx, k dz) \ .
		\end{equation*}
		and the 3-point operator for the 2nd derivatives in time leads to the extrapolation scheme
		\begin{block}{}
		\begin{equation*}
			 \frac{p_{j,k}^{n+1} - 2 p_{j,k}^n + p_{j,k}^{n-1}}{\mathrm{d}t^2} \ = \ c_j^2 ( \partial_x^2 p + \partial_z^2 p) \ + s_{j,k}^n
		\end{equation*}
		\end{block}
		where
		\begin{equation*}
			\begin{split}
				\partial_x^2 p \ &= \ \frac{p_{j+1,k}^{n} - 2 p_{j,k}^n + p_{j-1,k}^{n}}{\mathrm{d}x^2} \\ 
				\partial_z^2 p \ &= \ \frac{p_{j,k+1}^{n} - 2 p_{j,k}^n + p_{j,k-1}^{n}}{\mathrm{d}z^2}  
			\end{split}
		\end{equation*}
	\end{frame}
		
% ---------------------------------------------------------	

	\begin{frame}{Example}	
		\begin{itemize}
			\item[] Investigating the behaviour of the wavefield
			\item[] Simulating P-wave propagation in a reservoir scale model with maximum velocity $c_{max}=5 km/s$ and minimum velocity $c_{min}=3 km/s$
			\item[] Dominant frequency is 20 Hz $\longrightarrow$ dominant wavelength $\lambda_{dom}=c/f_{dom} = 150m$
			\item[] Simulating a spatial domain of $5\times 5 km$ and using a grid point distance $dx=10m$ resulting in 15 grid points per wavelength for the dominant frequency
		\end{itemize}
	\end{frame}
				
% ---------------------------------------------------------	

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

\subsection{Numerical anisotropy}

	\begin{frame}{Numerical anisotropy}
		Starting with the description of a plane harmonic wave propagating in 2D with wavenumber vector \textbf{k} pointing in the direction of propagation
		\begin{equation*}
			p(x,z,t) \ = \ e^{i (\boldsymbol{kx} - \omega t )} \ = \ e^{i (k_x x + k_z z - \omega t )} \  .
		\end{equation*} 
		and using the discretization of our 2D problem to obtain
		\begin{equation*}
			p^n_{j,k} \ = \  e^{i (k_x j \mathrm{d}x + k_z k \mathrm{d}x  - \omega n \mathrm{d}t )} \ .
		\end{equation*}
	\end{frame}		
				
% ---------------------------------------------------------	

	\begin{frame}{Numerical anisotropy}
		Injecting the formulation into the finite-difference approximation of the source-free 2D acoustic wave equation and following the same steps as done for the 1D numerical dispersion analysis leads to the following relation for the numerical phase velocity in 2D (assuming dx = dz)
		\begin{equation*}
			c^{num} ( k_x, k_z) \ = \ \frac{2}{\lambda \mathrm{d}t} \sin^{-1} \left[ \frac{\mathrm{d}t^2}{\mathrm{d}x^2} c^2 
    \left(\sin \left( \frac{k_x \mathrm{d}x}{2}\right) + \sin \left(\frac{k_z \mathrm{d}x}{2} \right) \right) \right]   \ . 
\end{equation*}
This relation can be analysed as a function of propagation direction noting that 
		\begin{equation*}
	\boldsymbol{k} \ = \ 
    \left( \begin{array}{c}
        k_x \\
        k_z 
    \end{array} \right)
    \ = \ 
    \left( \begin{array}{c}
        |{\bf k}| \cos{\alpha} \\
        |{\bf k}| \sin{\alpha}
    \end{array} \right)    
		\end{equation*}
	\end{frame}
				
% ---------------------------------------------------------	

	\begin{frame}{Numerical anisotropy}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_num_anis.png}
		\end{figure}
	\end{frame}		
	
% ---------------------------------------------------------	

	\begin{frame}{Fault zone waves - Landers, CA}
		\begin{figure}
			\includegraphics[width=0.30\textwidth]{\figpath/fig_fd_landers.png}
		\end{figure}
	\end{frame}		
				
% ---------------------------------------------------------	

\subsection{Choosing the Right Simulation Parameters}

	\begin{frame}{Choosing the Right Simulation Parameters}
		 \begin{itemize}
		 	\item[] The simulation of fault-zone trapped waves
			\item[] Noting: Scalar acoustic wave equation is mathematically identical to the SH-wave propagation problem
			\item[] Investigating the effects of a narrow 200m wide fault zone with a 25\% velocity decrease
		\end{itemize}
		\begin{table}
{
\begin{tabular}[t]{lll}
\ & Description & Value \\
\hline
$f_0$,$f_{max}$& Dominant, maximum frequency & 10 Hz, 30 Hz \\
$c_{min}$,$c_{max}$& Min., max. velocity & 2250 m/s, 3000 m/s \\
$x_{max}$,$z_{max}$& Min., max. extension & 10 km, 10 km \\
$t_{max}$ & Seismogram length & 3.5 s \\
\end{tabular}
}
\end{table}
	\end{frame}	
					
% ---------------------------------------------------------	

	\begin{frame}{Choosing the Right Simulation Parameters}
		\begin{itemize}[<+->]
			\item What is the smallest wavelength propagating in the medium? 
			\item[] {\it Answer:} With the $f_{max}=30 Hz$ and the smallest velocity $c_{min}=2250 m/s$ the shortest wavelength is given by $\lambda_{min}=c_{min}/f_{max}=75m$.
			\item How many numbers per smallest wavelength is required by the numerical method (given the wave propagation distance that needs to be covered)? 
			\item[] {\it Answer:} The dominant wavelength in the low-velocity medium is $\lambda_{dom}=c_{min}/f_{dom}=225m$. Thus we expect to propagate more than 20 wavelengths to the surface. As we use a 5-point operator we choose 20 points per dominant wavelength, resulting in about 7.5 points per smallest wavelength. 
		\end{itemize}
	\end{frame}	
					
% ---------------------------------------------------------	

	\begin{frame}{Choosing the Right Simulation Parameters}
		\begin{itemize}[<+->]
			\item What is the (in our case constant) grid spacing that needs to be implemented? 
			\item[] {\it Answer:} With 20 points per dominante wavelength we obtain $dx=\lambda_{dom}/20=11.25 m$ grid spacing. 
			\item What is the size of the physical domain and how many overal grid points are required?  
			\item[] {\it Answer:} The spatial extent of the 2D model is $10 km \times 10 km$. In each dimension we thus need $10000 m / dx \approx 900$ grid points leading to $900^2$ overall grid points.    
		\end{itemize}
	\end{frame}	
						
% ---------------------------------------------------------	

	\begin{frame}{Choosing the Right Simulation Parameters}
		\begin{itemize}[<+->]
			\item Will the seismogram(s) be influenced by artificial boundary reflections (i.e., is it necessary to implement absorbing boundaries or increase the model size)? 
			\item[] {\it Answer:} With the model setup as chosen - the fault zone in the middle of the model -  we will not expect problems with reflections from the boundaries. Thus no need to implement absorbing boundaries.  
			\item What is the maximum velocity in the model, and the resulting time step (given the grid increment and the CFL criterion)? 
			\item[] {\it Answer:} The maximum velocity in the model is $c_{max}=3000 m/s$. Assuming a CFL value of $\epsilon=0.7$ we can determine the time step required for a stable simulation as $dt = \epsilon dx/c_{max} = 0.0026 s$. 
		\end{itemize}
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Choosing the Right Simulation Parameters}
		\begin{itemize}[<+->]
			\item What is the overall number of time steps to be propagated?
			\item[] {\it Answer:} For a desired simulation time of $t_{max}=3.5 s$ the number of time steps required is $t_{max}/dt \approx 1300$. 
			\item How much core memory (RAM) will the simulation approximately require? 
			\item[] {\it Answer:} For a simple estimate we focus on the space-dependent fields that will constitute the largest part of the memory allocation. Those fields are 1) the velocity model, and 2) the pressure field at three different time levels,  and 3) two temporary fields containing the 2nd space derivatives. Assuming double precision floating point numbers (8 bytes per number) this will require approximately $6 \times 900^2 \time 8 bytes \approx 40 MBytes$.
		\end{itemize}
	\end{frame}	
							
% ---------------------------------------------------------	

	\begin{frame}{Choosing the Right Simulation Parameters}
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
				\includegraphics[width=0.8\textwidth]{\figpath/fig_ac2d_fz.png}
			\end{figure}		
		\end{minipage}
    \mbox{ }
    \hfill	
		\begin{minipage}{0.55\textwidth}
			\begin{itemize}
				\item[] Snapshot of pressure amplitude for a source located at the left edge of a fault zone model
				\item[] It indicates the occurence of head waves and the development of a dispersive wavefield trapped inside the low-velocity zone
			\end{itemize}
		\end{minipage}
	\end{frame}	
								
% ---------------------------------------------------------	

\section{Elastic Wave Propagation in 1D}	


% ---------------------------------------------------------
	
	\begin{frame}{Elastic Wave Propagation in 1D}
		In chapter 3 we introduced the stress-strain relation as
		\begin{equation*}
			\sigma_{ij} \ = \ \lambda \epsilon_{kk} \delta_{ij} + 2 \mu \epsilon_{ij}
		\end{equation*}
		where $\lambda, \mu$ are the Lam\'e parameters and $ \sigma_{ij}$ and  $\epsilon_{ij} $ are the stress and strain tensors \\
		Considering a 1D wave propagation in x-direction and particle motion in a horizontal direction y this relation reduces to
		\begin{equation*}
			\sigma_{xy} \ = \sigma_{yx} \ = \ 2 \mu \tfrac{1}{2} \left(  \partial_x u_y + \partial_y \underbrace{u_x}_{= 0} \right) \ = \ \mu \partial_x u_y   
		\end{equation*}
	\end{frame}
									
% ---------------------------------------------------------	

\subsection{Displacement - Stress Formulation}
	
	\begin{frame}{Displacement - Stress Formulation}
		\begin{block}{1D elastic wave equation}
			\begin{equation*}
				\rho \ddot{u} \ = \ \partial_x (\mu \, \partial_x) u + f
			\end{equation*}
		\end{block}
		\vspace{1cm}
		All space-dependent fields are defined at time level j through
		\begin{equation*}
			u_i^j \ = \ u(i \mathrm{d}x, j\mathrm{d}t)
		\end{equation*}
		\begin{figure}
			\begin{tikzpicture}[scale = .3]
\draw [fill] (0.2,0.2) circle [radius=0.1];
\node [below]  at (0.25,0) {$i-2$};
\draw [fill] (3.2,0.2) circle [radius=0.1];
\node [below]  at (3.25,0) {$i-1$};
\draw [fill] (6.2,0.2) circle [radius=0.1];
\node [below]  at (6.25,0) {$i$};
\draw [fill] (9.2,0.2) circle [radius=0.1];
\node [below]  at (9.25,0) {$i+1$};
\draw [fill] (12.2,0.2) circle [radius=0.1];
\node [below]  at (12.25,0) {$i+2$};
\end{tikzpicture}
		\end{figure}
	\end{frame}
									
% ---------------------------------------------------------	

	\begin{frame}{Displacement - Stress Formulation}
		Using the centred finite-difference method in space and time, the derivative of the displacement $u^\prime$ is
		\begin{equation*}
			u^\prime \ = \ \partial_x u \ = \ \frac{u_{i+1} - u_{i-1} }{2 \mathrm{d}x} 
		\end{equation*}
		To obtain the r.h.s. of 1D elastic wave equation we multiply by the shear modulus $\mu$
		\begin{equation*}
			\Eval{\mu \partial_x u}{i} \ = \Eval{\mu \ u^\prime}{i} \ = \ \mu_i \frac{u_{i+1} - u_{i-1} }{2 \mathrm{d}x}  
		\end{equation*}
	\end{frame}
										
% ---------------------------------------------------------	

	\begin{frame}{Displacement - Stress Formulation}		
		Taking the derivative of this equation results into
		\begin{equation*}
			\begin{split}
				\Eval{\partial_x \mu \partial_x u}{i} \ &= \frac{ \mu  u^\prime |_{i+1} - \mu  u^\prime |_{i-1}}{2 \mathrm{d}x} \\
  &= \ \frac{ \frac{\mu_{i+1} (u_{i+2} - u_i)}{2 \mathrm{d}x} - \frac{\mu_{i-1} (u_i - u_{i-2}) }{2 \mathrm{d}x}}
      {2   \mathrm{d}x} \\
  &= \ \frac{ \mu_{i+1}u_{i+2} - \mu_{i+1}u_i - \mu_{i-1}u_i + \mu_{i-1}u_{i-2} }{4 \mathrm{d}x^2}
			\end{split}
		\end{equation*}
		Approximating the l.h.s. of the wave equation with a centered scheme for the 2nd time derivative at time level $j$ leads to
		\begin{equation*}
			\begin{split}
				\rho_i \frac{ u_i^{j+1} - 2u_i^j + u_i^{j-1} }{\mathrm{d}t^2} \ &= \ \frac{ \mu_{i+1}u_{i+2}^j - \mu_{i+1}u_i^j - 
    \mu_{i-1}u_i^j + \mu_{i-1}u_{i-2}^j }{4 \mathrm{d}x^2} \\
    & \ \ \ + f_i^j
			\end{split}
		\end{equation*}
	\end{frame}
										
% ---------------------------------------------------------	

	\begin{frame}{Displacement - Stress Formulation}	
		Final extrapolation scheme for the displacement-stress 1D elastic wave equation using a central difference scheme
		\begin{equation*}
			\begin{split}
				u_i^{j+1} \ &= \ \frac{\mathrm{d}t^2}{4 \rho_i \mathrm{d}x^2} \left[ \mu_{i+1}u_{i+2}^j - \mu_{i+1}u_i^j - 
    \mu_{i-1}u_i^j + \mu_{i-1}u_{i-2}^j \right] \\
    &\ \ \ \ \ + 2u_i^j -u_i^{j-1} + \frac{\mathrm{d}t^2 }{\rho_i} f_i^j  
			\end{split}
		\end{equation*}
	\end{frame}
											
% ---------------------------------------------------------	

\subsection{Velocity - Stress Formulation}

	\begin{frame}{Velocity - Stress Formulation}
		Defining velocity v and stress component $\sigma$ as
		\begin{equation*}
			\begin{split}
				\partial_t u \ &= \ v \\
				\sigma \ &= \ \mu \partial_x u
			\end{split}
		\end{equation*}
		and assuming space-time dependencies leads to the wave equation
		\begin{equation*}
			\begin{split}
				\rho \partial_t v \ &= \ \partial_x \sigma + f \\
				\dot{\sigma} \ &= \ \mu \partial_x v  
			\end{split}
		\end{equation*}
		Our unknowns are
		\begin{equation*}
			v_i^j \ = \ v(i \mathrm{d}x, j\mathrm{d}t)
		\end{equation*}
	\end{frame}
																
%% ---------------------------------------------------------	
%
	%\begin{frame}{Velocity - Stress Formulation}
		%\centering Spatio-temporal grid staggering in 1D
	%
		%
		%\begin{figure}
			%\begin{tikzpicture}[scale = .7]
%\draw [step = 1, gray, very thin] (0,0) grid (6,6);
%\draw [->] (0, 0) -- (0,6.5);
%\node [left]  at (0,2) {$j-1$};
%\node [left]  at (0,3) {$j - \tfrac{1}{2}$};
%\node [left]  at (0,4) {$j$};
%\node [left]  at (0,5) {$j + \tfrac{1}{2}$};
%
%\draw [->] (0, 0) -- (6.5,0);
%\node [below]  at (2,0) {$i$};
%\node [below]  at (3,0) {$i+\tfrac{1}{2}$};
%\node [below]  at (4,0) {$i+1$};
%
%\draw (1,2) circle [radius = 0.05];
%\draw (1,4) circle [radius = 0.05];
%
%
%\draw [fill] (2,1) circle [radius = 0.05];
%\draw [fill] (2,3) circle [radius = 0.05];
%\draw [fill] (2,5) circle [radius = 0.05];
%
%\draw (3,2) circle [radius = 0.05];
%\draw (3,4) circle [radius = 0.05];
%
%\draw [fill] (4,1) circle [radius = 0.05];
%\draw [fill] (4,3) circle [radius = 0.05];
%\draw [fill] (4,5) circle [radius = 0.05];
%
%\draw (5,2) circle [radius = 0.05];
%\draw (5,4) circle [radius = 0.05];
%
%\draw [fill] (6,1) circle [radius = 0.05];
%\draw [fill] (6,3) circle [radius = 0.05];
%\draw [fill] (6,5) circle [radius = 0.05];
%
%\draw (1.5, 2.5) rectangle (3.5, 4.5);
%\node at (2.5, 3.5) {$i,j$};
%
%\draw [fill] (7,4.2) circle [radius = 0.05];
%\node [right] at (7,4.2) {$\, v, \rho$};
%\draw        (7,3.8) circle [radius = 0.05];
%\node [right] at (7,3.8) {$\, \mu , \sigma$};
			%\end{tikzpicture}
		%\end{figure}	
	%\end{frame}
	%
%% ---------------------------------------------------------	
%
	%\begin{frame}{Velocity - Stress Formulation}
		%The following computational scheme does the trick:
		%\begin{equation*}
			%\begin{split}
				%\frac{v_i^{j+ \tfrac{1}{2}} - v_i^{j- \tfrac{1}{2}} }{\mathrm{d}t} \ &= \ \frac{1}{\rho_i}
    			%\frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{\mathrm{d}x} + \frac{f_i^j}{\rho_i} \\
				%\frac{\sigma_{i+\tfrac{1}{2}}^{j+1} - \sigma_{i+\tfrac{1}{2}}^j }{\mathrm{d}t} 
    %\ &= \ \mu_{i+\tfrac{1}{2}}  \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }
    %{\mathrm{d}x}
			%\end{split}
		%\end{equation*}
		%leading to the extrapolation scheme
		%\begin{equation*}
			%\begin{split}
				%v_i^{j+ \tfrac{1}{2}} \ &= \ \frac{dt}{\rho_i}
    %\frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{\mathrm{d}x} \ + \ v_i^{j- \tfrac{1}{2}}  + \frac{dt}{\rho_i} \ f_i^j  \\
%\sigma_{i+\tfrac{1}{2}}^{j+1} \ &= \ \mathrm{d}t \ \mu_{i+\tfrac{1}{2}} \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }{\mathrm{d}x} \ + \ \sigma_{i+\tfrac{1}{2}}^j  
   		 	%\end{split}
		%\end{equation*}
	%\end{frame}
	
	
	% attempt tp merge

\begin{frame}{Velocity - Stress: Staggered grid formulation}
	
	\begin{minipage}{0.48\textwidth}	
		\begin{figure}
			\begin{tikzpicture}[scale = .7]
\draw [step = 1, gray, very thin] (0,0) grid (6,6);
\draw [->] (0, 0) -- (0,6.5);
\node [left]  at (0,2) {$j-1$};
\node [left]  at (0,3) {$j - \tfrac{1}{2}$};
\node [left]  at (0,4) {$j$};
\node [left]  at (0,5) {$j + \tfrac{1}{2}$};

\draw [->] (0, 0) -- (6.5,0);
\node [below]  at (2,0) {$i$};
\node [below]  at (3,0) {$i+\tfrac{1}{2}$};
\node [below]  at (4,0) {$i+1$};

\draw (1,2) circle [radius = 0.05];
\draw (1,4) circle [radius = 0.05];


\draw [fill] (2,1) circle [radius = 0.05];
\draw [fill] (2,3) circle [radius = 0.05];
\draw [fill] (2,5) circle [radius = 0.05];

\draw (3,2) circle [radius = 0.05];
\draw (3,4) circle [radius = 0.05];

\draw [fill] (4,1) circle [radius = 0.05];
\draw [fill] (4,3) circle [radius = 0.05];
\draw [fill] (4,5) circle [radius = 0.05];

\draw (5,2) circle [radius = 0.05];
\draw (5,4) circle [radius = 0.05];

\draw [fill] (6,1) circle [radius = 0.05];
\draw [fill] (6,3) circle [radius = 0.05];
\draw [fill] (6,5) circle [radius = 0.05];

\draw (1.5, 2.5) rectangle (3.5, 4.5);
\node at (2.5, 3.5) {$i,j$};

\draw [fill] (7,4.2) circle [radius = 0.05];
\node [right] at (7,4.2) {$\, v, \rho$};
\draw        (7,3.8) circle [radius = 0.05];
\node [right] at (7,3.8) {$\, \mu , \sigma$};
			\end{tikzpicture}
		\end{figure}	
		
		\end{minipage}
		\mbox{}
			\begin{minipage}{0.5\textwidth}	
		\begin{equation*}
			\begin{split}
				\frac{v_i^{j+ \tfrac{1}{2}} - v_i^{j- \tfrac{1}{2}} }{\mathrm{d}t} \ &= \ \frac{1}{\rho_i}
    			\frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{\mathrm{d}x} + \frac{f_i^j}{\rho_i} \\
				\frac{\sigma_{i+\tfrac{1}{2}}^{j+1} - \sigma_{i+\tfrac{1}{2}}^j }{\mathrm{d}t} 
    \ &= \ \mu_{i+\tfrac{1}{2}}  \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }
    {\mathrm{d}x}
			\end{split}
		\end{equation*}
		leading to the extrapolation scheme
		\begin{equation*}
			\begin{split}
				v_i^{j+ \tfrac{1}{2}} \ &= \ \frac{dt}{\rho_i}
    \frac{\sigma_{i + \tfrac{1}{2}}^j - \sigma_{i - \tfrac{1}{2}}^j }{\mathrm{d}x} \ + \ v_i^{j- \tfrac{1}{2}}  + \frac{dt}{\rho_i} \ f_i^j  \\
\sigma_{i+\tfrac{1}{2}}^{j+1} \ &= \ \mathrm{d}t \ \mu_{i+\tfrac{1}{2}} \frac{v_{i + 1}^{j +\tfrac{1}{2}} - v_i^{j + \tfrac{1}{2}} }{\mathrm{d}x} \ + \ \sigma_{i+\tfrac{1}{2}}^j  
   		 	\end{split}
		\end{equation*}
		
		\end{minipage}
	
		
	\end{frame}
	
	
% ---------------------------------------------------------	

	\begin{frame}{Example}
		\begin{minipage}{0.55\textwidth}
			Simulation parameters and Matlab code fragment for 1D velocity-stress simulation \\
			\begin{table}
			{
			\begin{tabular}[t]{ll}
				Parameter & Value \\
				\hline
				$\rho$ & 2500 $kg/m^3$ \\
				$\mu $ & 5 $\times$ 50 GPa \\
				$ v_S$ & 4500 $m/s$ \\
				$x_{max}$ & $1000$ km \\
				dx & 1000 m \\
				dt & 0.18 s \\
				$f_0$ & 1/15 Hz \\ 
			\end{tabular}
			}
			\end{table}
			\vspace{2cm}
		\end{minipage}
	 \mbox{ }
 		\begin{minipage}{0.3\textwidth}
	\begin{figure}
			\includegraphics[width=1.85\textwidth]{python_staggered.png}
		\end{figure}
		\end{minipage}
	\end{frame}
				
% ---------------------------------------------------------	

	\begin{frame}{Example}
		\begin{figure}
			\includegraphics[width=0.5\textwidth]{\figpath/fig_el1d.png}
		\end{figure}
	\end{frame}				
					
% ---------------------------------------------------------	
	
\subsection{Velocity - Stress: Dispersion}

	\begin{frame}{Velocity - Stress: Dispersion}
	The condition for stable calculations can be obtained by
		\begin{equation*}
			\sin \left( \frac{\omega \mathrm{d}t}{2} \right) \ = \ \pm \sqrt{ \frac{\mu}{\rho}}
    		\, \frac{\mathrm{d}t}{\mathrm{d}x} \sin \left( \frac{k \mathrm{d}x}{2} \right) 
		\end{equation*}
		For the numerical phase velocity we obtain
		\begin{equation*}
			c^{num} \ = \ \frac{\omega}{k} \ = \ \frac{\lambda}{\pi \mathrm{d}t} 
   			 \sin^{-1} \left( c_0 \frac{\mathrm{d}t}{\mathrm{d}x} \sin{\frac{\pi \mathrm{d}x}{\lambda}} \right)
		\end{equation*}
		where $k=2 \pi/\lambda$ \\
		The group velocity is given as
		\begin{equation*}
			\frac{\mathrm{d}\omega}{\mathrm{d}k} \ = \ \frac{c_0 \cos{\frac{\pi \mathrm{d}x}{\lambda}} }
    { \left[ 1 - \left( c_0 \frac{\mathrm{d}t}{\mathrm{d}x} 
    \sin{ \frac{\pi \mathrm{d}x}{\lambda} } \right)^2 \right]^{\tfrac{1}{2}}}
		\end{equation*}
	\end{frame}
					
% ---------------------------------------------------------	

	\begin{frame}{Velocity - Stress: Dispersion}
		\begin{minipage}{0.45\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_el1d_courant1.png}
		\end{figure}		
		\end{minipage}
	 \mbox{ }
    \hfill		
		\begin{minipage}{0.45\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_el1d_courant2.png}
		\end{figure}		
		\end{minipage}
	\end{frame}
					
% ---------------------------------------------------------	

\subsection{Elastic 2D - Staggered Grids}

	\begin{frame}{Elastic 2D - Staggered Grids}
		\begin{minipage}{0.55\textwidth}
			\begin{figure}
				\begin{tikzpicture}[scale = .5]
\draw [step = 1, gray, very thin] (0,0) grid (6,6);
\draw [->] (0, 0) -- (0,6.5);
\node [left]  at (0,6) {$z$};
\draw [->] (0, 0) -- (6.5,0);
\node [below]  at (6,0) {$x$};

\draw (1,1) circle [radius = 0.08];
\draw (0.93,1.93) rectangle (1.08, 2.08);
\draw (1,3) circle [radius = 0.08];
\draw (0.93,3.93) rectangle (1.08, 4.08);
\draw (1,5) circle [radius = 0.08];

\draw [fill] (2,1) circle [radius = 0.08];
\draw (1.9,1.9)--(2.1,1.9)--(2, 2.1)--(1.9,1.9) ; % triangle
\draw [fill] (2,3) circle [radius = 0.08];
\draw (1.9,3.9)--(2.1,3.9)--(2, 4.1)--(1.9,3.9) ; % triangle
\draw [fill] (2,5) circle [radius = 0.08];

\draw (3,1) circle [radius = 0.08];
\draw (2.93,1.93) rectangle (3.08, 2.08);
\draw (3,3) circle [radius = 0.08];
\draw (2.93,3.93) rectangle (3.08, 4.08);
\draw (3,5) circle [radius = 0.08];

\draw [fill] (4,1) circle [radius = 0.08];
\draw (3.9,1.9)--(4.1,1.9)--(4, 2.1)--(3.9,1.9) ; % triangle
\draw [fill] (4,3) circle [radius = 0.08];
\draw (3.9,3.9)--(4.1,3.9)--(4, 4.1)--(3.9,3.9) ; % triangle
\draw [fill] (4,5) circle [radius = 0.08];

\draw (5,1) circle [radius = 0.08];
\draw (4.93,1.93) rectangle (5.08, 2.08);
\draw (5,3) circle [radius = 0.08];
\draw (4.93,3.93) rectangle (5.08, 4.08);
\draw (5,5) circle [radius = 0.08];

\draw (1.5,2.5) rectangle (3.5, 4.5);
\node at (2.5,3.5) {$i,j$};

\draw [fill] (6.5,5.5) circle [radius = 0.08];
\node [right] at (6.5,5.5) {$\, v_z, \rho$};
\draw (6.43,5.03) rectangle (6.58, 5.18);
\node [right] at (6.5,5.1) {$\, v_x, \rho$};
\draw (6.4,4.6)--(6.6,4.6)--(6.5, 4.8)--(6.4,4.6) ; % triangle
\node [right] at (6.5,4.7) {$\, \dot{\sigma}_{zz}, \dot{\sigma}_{xx}, \lambda, \mu$};
\draw (6.5,4.3) circle [radius = 0.08];
\node [right] at (6.5,4.3) {$\, \dot{\sigma}_{xz}, \mu$};
				\end{tikzpicture}
			\end{figure}		
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{enumerate}
					\item Contains fundamental aspects of grid staggering for the stress-strain relation for higher dimensions
					\item The extension to 3D is straight forward
					\item Scheme in seismology most widely used
				\end{enumerate}
			\end{block}
		\end{minipage}
	\end{frame}
					
% ---------------------------------------------------------	

	\begin{frame}{Elastic 2D - Staggered Grids}
		The stress-strain relation in 1D
		\begin{equation*}
			\dot{\sigma}_{ij} \ = \ \lambda \dot{\epsilon}_{kk} \delta_{ij} + 2\mu \dot{\epsilon}_{ij} 
		\end{equation*}
		and rewriting it in 2D using the definition of the strain tensor to obtain for each component
		\begin{equation*}
			\begin{split}
				\dot{\sigma}_{xx} \ &= \ ( \lambda + 2\mu ) \partial_x v_x + \lambda \partial_z v_z \\
				\dot{\sigma}_{zz} \ &= \ ( \lambda + 2\mu ) \partial_z v_z + \lambda \partial_x v_x \\
				\dot{\sigma}_{xz} \ &= \ \mu ( \partial_x v_z + \lambda \partial_z v_x )
			\end{split}
		\end{equation*}
	\end{frame}
	
	
% ---------------------------------------------------------	

\subsection{Elastic 3D - Staggered Grids}

	\begin{frame}{Elastic 3D - Staggered Grids}
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_3D.png}
		\end{figure}	
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{enumerate}
					\item Most common FD scheme today
					\item Usually 4-point operators
					\item Time extrapolation with Runge-Kutta
				\end{enumerate}
			\end{block}
		\end{minipage}
	\end{frame}


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

	
\subsection{Heterogeneous Earth models}

	\begin{frame}{Layered models}
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_interface.png}
		\end{figure}	
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item Layers might not coincide with grid lines 
					\item Blocky representation of interfaces (errors)
					\item Grid stretching, curvi-linear grids, homogenization
				\end{itemize}
			\end{block}
		\end{minipage}
	\end{frame}


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

	\subsection{Optimizing Operators}

	\begin{frame}{Optimal Operators (Geller et al.)}
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_opt.png}
		\end{figure}	
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item FD operators spread in time and space
					\item Implicit scheme turned explicit using Born approximation
				\end{itemize}
			\end{block}
		\end{minipage}
	\end{frame}

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

	\subsection{Optimizing Operators}

	\begin{frame}{Optimal Operators (Geller et al.)}
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_geller.png}
		\end{figure}	
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{block}{}
				\begin{itemize}
					\item Substantial accuracy improvement
					\item Negligible extra costs
				\end{itemize}
			\end{block}
		\end{minipage}
	\end{frame}			
			
% ---------------------------------------------------------	

	\subsection{Unstructured Grids}
\begin{frame}{Unstructured grids}
\begin{figure}
			\includegraphics[width=0.65\textwidth]{\figpath/fig_fd_voronoi1.png}
		\end{figure} 
		
	\end{frame}	
	
	
% ---------------------------------------------------------	

	\subsection{Other coordinate systems}
\begin{frame}{Spherical coordinates}
		\begin{minipage}{0.55\textwidth}
		\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_axi4.png}
		\end{figure}	
		\end{minipage}
	 \mbox{ }
    \hfill				
		\begin{minipage}{0.4\textwidth}
			\begin{figure}
			\includegraphics[width=0.95\textwidth]{\figpath/fig_fd_shaxi.png}
		\end{figure}	
		\end{minipage}
		
	\end{frame}	
					
% ---------------------------------------------------------	

	\begin{frame}{Summary}
	\scriptsize
	\begin{itemize}
		\item{Replacing the partial derivatives by finite differences allows partial differential equations such as the wave equation to be solved directly for (in principle) arbitrarily heterogeneous media.}
		\item{The resulting space-time discretization leads to unphysical phenomena such as numerical dispersion that can only be avoided by sampling with enough number of grid points per wavelength.}
		\item{The accuracy of finite-difference operators can be improved by using information from more grid points (i.e., longer operators). The weights for the grid points can be obtained using Taylor series.}  
		\item{Classic plane-wave analysis of the approximative scheme leads to the famous Courant-Friedrich-Levy (CFL) criterion that restricts the choice of the space-time discretization.}
		\item{In higher dimensions the error of the wave propagation becomes anisotropic. In regular-spaced grids the most accurate direction is at $45^\circ$ w.r.t. the grid axes.}
		\item{The implementation of boundary conditions in the case of finite-differences needs special care.}	
	\end{itemize}
	\end{frame}
					
% ---------------------------------------------------------	

	\begin{frame}{Comprehensive Questions}
	\scriptsize
		\begin{enumerate}
			\item[1] Are finite-difference (FD) based approximations of partial-differential equations unique (give arguments)?
			\item[2] What strategies are there to improve the accuracy of finite-difference based derivatives? Give the procedures in words.
			\item[3]	What is stability in connection with FD algorithms, give the relevant condition for the 1D wave propagation problem. 
			\item[4] What is convergence? 
			\item[5] What is the difference between physical and numerical dispersion?
			\item[6] Which propagation direction is most accurate on a rectangular (square) grid, any reasons?
			\item[7] How would you check whether an FD simulation is accurate (a - homogeneous medium, b - strongly heterogeneous medium)? 
		\end{enumerate}
	\end{frame}
					
% ---------------------------------------------------------	

	\begin{frame}{Theoretical problems}
	\scriptsize
		\begin{enumerate}
\item[8]  
Show that the  5-point operator given in Eq.~\ref{5-pt-fd} is an approximation for the second derivative of $f(x)$ with respect to $x$ at position $x$. Hint: Use Taylor's theorem. What is the leading order of the error term?
\item[9]
A seismometer consists of a spring with damping parameter $\epsilon$, and eigenfrequency $\omega_0$. The seismometer is excited by the (given) ground motion $\ddot{u}(t)$. The relative motion of the seismometer mass x(t) is governed by the following equation
\begin{displaymath}
\ddot{x}+2\epsilon \dot{x} + \omega^2_0 x = \ddot{u}
\end{displaymath}
Replace the derivatives on the l.h.s. with finite differences. Solve for $x(t+dt)$. Note: A good strategy in this example is to center the differences at the same point in time.  
			\end{enumerate}
	\end{frame}
			
% ---------------------------------------------------------	

	\begin{frame}{Theoretical problems}
	\scriptsize
		\begin{enumerate}
\item[10] 
Certain isotopes (e.g., $_9$Be) are washed into the sea by rivers and then mixed by advection through ocean currents and diffusion. In addition, the isotopes are removed from the system through biomechanical processes (e.g., death). These processes can be described through the diffusion-advection-reaction equation (concentration $C(x,t)$, diffusivity $k$ (const), reactivity $R(x)$, source $p(x)$, velocity $v(x)$). Substitute in the 1-D equation below the partial differentials with finite differences and extrapolate to $C(t+dt)$: 
\begin{displaymath}
\partial_t C = k \partial^2_x C + v_x \partial_x C - R C + p
\end{displaymath}
How could a "ring-current" be simulated with this 1-D equation mimicking an oceanic gyre?   
			\end{enumerate}
	\end{frame}
			
% ---------------------------------------------------------	

	\begin{frame}{Theoretical problems}
	\scriptsize
		\begin{enumerate}
\item[11]
You want to simulate 2-D acoustic wave propagation in a medium with size 1000km x 1000km. You want to model wave propagation up to a period of 10 seconds. The maximum velocity c is 8km/s, the minimum velocity is 4km/s. Your numerical algorithm requires 20 grid points per wavelength to be accurate for the propagation distances of interest. What space increment dx do you need for the simulation?
The stability criterion says that maximum velocity c, space increment dx and time increment dt are related by $const=c dt/dx$. You want a seismogram length of 500 seconds. How many time steps do you have to simulate, when const=0.5? 
\end{enumerate}
	\end{frame}

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

\end{document}