diff git a/AUTHORS b/AUTHORS
index 1743944e44f22e3f818f865b69661450243f4939..fc92491a1181a8ef39ec0dd5c6fa9b55715f440c 100644
 a/AUTHORS
+++ b/AUTHORS
@@ 7,4 +7,5 @@ Olaf Hellwig
Stefan Jetschny
Daniel Koehn
Andre Kurzmann
Martin Schaefer
\ No newline at end of file
+Martin Schaefer
+Florian Wittkamp
\ No newline at end of file
diff git a/doc/guide_sofi3D/tex/eps/error_1D_courant.pdf b/doc/guide_sofi3D/tex/eps/error_1D_courant.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..577b1fb789aa782665acaf7f7412530ee40077b2
Binary files /dev/null and b/doc/guide_sofi3D/tex/eps/error_1D_courant.pdf differ
diff git a/doc/guide_sofi3D/tex/eps/error_3D_courant.pdf b/doc/guide_sofi3D/tex/eps/error_3D_courant.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..565ca47a577ca19d3f66e6cbd349c03b3bc3f747
Binary files /dev/null and b/doc/guide_sofi3D/tex/eps/error_3D_courant.pdf differ
diff git a/doc/guide_sofi3D/tex/eps/seismo.pdf b/doc/guide_sofi3D/tex/eps/seismo.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..caa15eb11742fe8751b1388be8397fd446d577f5
Binary files /dev/null and b/doc/guide_sofi3D/tex/eps/seismo.pdf differ
diff git a/doc/guide_sofi3D/tex/guide_sofi3D.tex b/doc/guide_sofi3D/tex/guide_sofi3D.tex
index 402212b467a167c9365399ad71156b6bdd004c0b..d297915d5f39b9f2aeb182b5a3aa07fd423deb42 100644
 a/doc/guide_sofi3D/tex/guide_sofi3D.tex
+++ b/doc/guide_sofi3D/tex/guide_sofi3D.tex
@@ 323,13 +323,127 @@ The resulting coefficients are $\beta_1=9/8$ and $\beta_2=1/24$, so the forward
The coefficients $\beta_i$ in the FD operator are called {\bf{Taylor coefficients}}. The accuracy of higher order FD operators can be improved significantly simply by slightly changes in the FD coefficients \cite{holberg:87}. These numerically optimized coefficients are called {\bf{Holberg coefficients}}. An accuracy comparison between Taylor and Holberg coefficients can be found in section \ref{comp_taylor_holberg}.
\subsection{Higher order temporal FDoperators}
To achieve a better performance and higher accuracy of the FD simulation higher order temporal FDoperators can be chosen. The order of the temporal integration is increased by the AdamsBashforth method (AB).
The AB method uses previous time levels of the wave field which will be weighted by coefficients \cite{ghrist2000staggered}.\\ \\
Due to previous time levels have to be saved, memory usage will by increase by a factor of 2.3 (elastic) and 2.6 (elastic) for third order and fourth order, respectively. In addition the calculation time will be slightly higher for higher temporal orders.\\\\
However the numerical dispersion of the AB method is much lower than that of the widely used secondorder leapfrog method.
Numerical dissipation is introduced by the AB method which is significantly smaller for the AB method of fourthorder accuracy then of thirdorder accuracy. Moreover the stability limit will be lower for the higher order temporal orders (see \ref{courandt}).\\ \\
A detailed analysis of the AdamsBashforth method applied to seismic wave simulation with FDTD is given by \cite{bohlen2015higher}.
+\subsection{AdamsBashforth higherorder accurate time integrators}
+To improve the precision of the temporal discretization the AdamsBashforth thirdorder and fourthorder accurate time integrator can be used.
+The staggered AdamsBashforth method (ABS) is a multistep method that uses previously calculated wavefields to increase the accuracy order in time.\\
+In \cite{bohlen2015higher2} and \cite{bohlen2015higher} we give a detailed study of the performance of the ABS method. The analysis shows that the numerical dispersion is much lower than that of the widely used secondorder leapfrog method. Numerical dissipation is introduced by the ABS method which is significantly smaller for the ABS method of fourthorder accuracy. In 1D and 3D simulation experiments the convincing improvements of simulation accuracy of the fourthorder ABS method is verified. In a realistic elastic 3D scenario the computing time reduces by a factor of approximately 2.2, whereas the memory requirements increase by approximately the same factor. \\
+The ABS method thus provides an alternative strategy to increase the time accuracy by investing computer memory instead of computing time.\\
+
+In the following we provide an extract of \cite{bohlen2015higher2} to demonstrate the underlying theory. Additionally seismograms of different temporal orders of accuracy and an analysis of the numerical simulation error are shown.\\
+\subsubsection{Theory}
+For the sake of simplicity we illustrate the staggered AdamsBashforth method (ABSmethod) using the 1D acoustic wave equation in velocitystress formulation
+\begin{eqnarray}
+\frac{\partial p(x,t)}{\partial t}&=&  \pi (x) \frac{\partial v(x,t)}{\partial x} \notag \\
+\frac{\partial v(x,t)}{\partial t}&=& \rho^{1}(x) \frac{\partial p(x,t)}{\partial x}
+\label{waveequation}
+\end{eqnarray}
+The wavefield variables are the pressure $p(x,t)$ and the particle velocity $v(x,t)$. The material is described by the Pwave modulus $\pi (x)$ and the mass density $\rho(x)$. For simplicity we omit the temporal ($t$) and spatial dependencies ($x$) in the following.
+
+Using the conventional second order staggered grid approximation to the first order time derivative we obtain
+\begin{eqnarray}
+\frac{p ^{n+1/2}p^{n1/2}}{\Delta t}&=&  \pi \frac{\partial v}{\partial x}\bigg ^n + \mathcal{O}(\Delta t^2) \nonumber \\
+\frac{v^nv^{n1}}{\Delta t}&=& \rho^{1} \frac{\partial p}{\partial x}\bigg ^{(n1/2)} + \mathcal{O}(\Delta t^2)
+\label{2ndorder}
+\end{eqnarray}
+
+This results in the conventional explicit second order accurate timestepping (leapfrog) scheme
+\begin{eqnarray}
+p^{n+1/2}&=&p^{n1/2}  \Delta t \pi \frac{\partial v}{\partial x}\bigg ^n + \mathcal{O}(\Delta t^2) \nonumber \\
+v^n&=&v^{n1} \Delta t \rho^{1} \frac{\partial p}{\partial x}\bigg ^{(n1/2)} + \mathcal{O}(\Delta t^2)
+\label{2ndorderexplicit}
+\end{eqnarray}
+which requires no additional storage of the wavefield variables $p$ and $v$.
+
+With the ABSmethod the order of the temporal integration can be increased by using previous time levels of the right hand sides of equation \ref{2ndorder} \cite{ghrist2000staggered}.
+The ABSmethod thus requires the storage of previous time levels of spatial derivatives of the pressure $p$ and particle velocity $v$. The time update for the ABSmethod thus reads
+\begin{eqnarray}
+p^{n+1/2}&=&p^{n1/2}  \Delta t \pi \sum_{k=0}^{M1} a_k \frac{\partial v}{\partial x}\bigg ^{nk} + \mathcal{O}(\Delta t^M) \nonumber \\
+v^n&=&v^{n1} \Delta t \rho^{1} \sum_{k=0}^{M1} a_k \frac{\partial p}{\partial x}\bigg ^{(n1/2k)} + \mathcal{O}(\Delta t^M)
+\label{absmethod}
+\end{eqnarray}
+The weights for the time accuracy orders $M=2,3,4$ are given in Table \ref{absweights}.
+
+\begin{table}[ht]
+\begin{center}
+\begin{tabular}{ccccc}\hline
+$M$ & $a_0$ & $a_1$ & $a_2$ & $a_3$ \\ \hline
+2 & 1 & 0 & 0 & 0 \\
+3 & 25/24 & 1/12 & 1/24 & 0 \\
+4 & 13/12 & 5/24 & 1/6 & 1/24 \\ \hline
+\end{tabular}
+\end{center}
+\caption{\label{absweights}Weights used for the summation of previous time levels in the AdamsBashforth method \protect\cite{ghrist2000staggered}.}
+\end{table}
+\newpage
+\subsubsection{Accuracy in 1D}
+We first compare seismograms calculated with the 1D ABSmethod using equations \ref{absmethod} for different orders of accuracy in time ($M$). We use a 1D homogeneous acoustic
+medium with a wave velocity of $c=3500$\,m/s and constant density of $\rho=2000$\,kg/m$^3$ . The source signal is a Ricker signal with a center frequency of 600\,Hz. We choose a large sourcereceiver distance of 120 dominant wavelength to emphasize the effects of numerical dispersion and numerical dissipation in the synthetic seismograms. The spatial derivatives are
+computed with high accuracy using a centered staggered FD stencil of 8th order accuracy. The spatial grid spacing is held constant at $\Delta x=0.4$\,m corresponding to approximately 14 grid points per dominant wavelength. Discrepancies to the analytical solution (timeshifted Ricker signal) are thus mainly caused by the chosen time step interval $\Delta t$ and the order of the temporal discretization ($M=2,3,4$). The numerical results for Courant numbers $r=c \Delta t/\Delta x=0.4, 0.2, 0.1$ and temporal orders of accuracy $M=2, 3, 4$ are compared with the analytical solution in Figure \ref{fig_1dfdseis}. For a large Courant number of $r=0.4$ (corresponding to a large time step interval) and second order approximation of the time derivative ($M=2$) we can observe a large time dispersion error causing a leading phase with high amplitude (Figure \ref{fig_1dfdseis}, topleft). Figure \ref{fig_1dfdseis} compares two ways to reduce this time discretization error. The conventional way is to stay with the given temporal accuracy order (typically $M=2$) and then reduce the time step interval, i.e. Courant number. The resulting waveforms are shown in the rows of Figure \ref{fig_1dfdseis}.
+Reducing the Courant number (time step interval) will increase the computation time proportional to $1/r$. The alternative way proposed in this study is to increase the order of the temporal discretization with the ABSmethod which requires to store $M1$ previously calculated spatial derivates of wavefields. The resulting waveforms are shown in the columns of Figure \ref{fig_1dfdseis}. \\ \\
+%
+%In Figure \ref{fig_1dfdseis} we see that the simulation accuracy increases with decreasing Courant number and increasing accuracy order. In order to quantify the corresponding numerical simulation error we calculated the normalized L2misfit between the numerical and analytical solution for different Courant numbers and accuracy orders
+%$M = 2, 3, 4$.
+%In Figure \ref{fig_L2} we plot the normalized L2misfit over the number of time steps $NT=T/\Delta t=T c/(r \Delta x)$ where the propagation time of waves is held fixed at $T=0.24$\,s. We can observe the expected higher order reduction of the simulation accuracy which reduces proportional to $\Delta t^{M}$ or $NT^{M}$. For a given accuracy level of $E=0.1\%$ (horizontal line in Figure \ref{fig_L2}) the required numbers of time steps are 39233, 13938, and 8704 for $M=2,M=3$ and $M=4$, respectively. The number of time steps thus reduces to $36\%$ and $22\%$ when we increase the temporal order from $M=2$ to $M=3$ and $M=4$, respectively (Figure \ref{fig_L2} ). The ABS methods thus allows to significantly decrease the number of time steps and thus to save computing time. \\ \\
+\begin{figure*}[!h]
+\begin{center}
+\includegraphics[width=\textwidth]{eps/seismo.pdf}
+\caption{Seismograms (red lines) calculated for the 1D case for different Courant numbers $r$ and accuracy orders $M$. $M=2$ corresponds to the classical second order leapfrog scheme (equations \ref{2ndorder}). $M=3,4$ correspond to the multistep ABS method (equations \ref{absmethod}). The analytical solution is plotted as a black line. }
+\label{fig_1dfdseis}
+\end{center}
+\end{figure*}
+%\begin{figure}[h]
+%\begin{center}
+%\includegraphics[scale=0.7]{eps/error_1D_courant}
+%\caption{Relative error of 1D simulations (Figure \ref{fig_1dfdseis}).
+%The normalized L2 norm is plotted over the required number of time steps ($NT$). The temporal accuracy orders are $M=2,3,4$. The spatial accuracy is fixed at order $N=8$.
+%The number of time steps to achieve a given accuracy of $E=0.01\%$ (horizontal line) reduce to $36\%$ and $22\%$ for the temporal orders $M=3$ and $M=4$, respectively.
+%For large $NT$ (small time step intervals) the error converges to the fixed error of the spatial discretization. }
+%\label{fig_L2}
+%\end{center}
+%\end{figure}
+
+\newpage
+\subsubsection{Accurace in 3D}
+In Figure \ref{fig_1dfdseis} we see that the simulation accuracy in 1D increases with decreasing Courant number and increasing accuracy order. In order to quantify the corresponding numerical simulation error of the 3D elastic update scheme we calculated the normalized L2misfit between the numerical and analytical solution for different Courant numbers and accuracy orders
+$M = 2, 3, 4$.
+We compared the synthetic seismograms with an analytical solution for an explosive point source in a homogeneous full space. The elastic model parameters are $v_p=3500$\,m/s
+and $v_s=2000$\,m/s for the P and Svelocities, respectively, and $\varrho=2000$\,kg/m$^3$ for the density. The explosive point source generates a Ricker signal with a dominant frequency of $600$\,Hz. We analyze the accuracy of the direct Pwave in a distance of $180$\,m corresponding to 30 dominant wavelength. The size of the model grid is 800x400x400 grid points. The direct Pwave is spatially sampled with approximately 14 grid points per dominant wavelength. The influence of numerical dispersion due to the discretization in space is small because of the
+chosen high spatial accuracy order of $N=8$. The simulations are performed on 80 cores on a smallscale shared memory cluster.
+
+As a measure of accuracy we use the normalized L2norm between the numerical and analytical seismograms which only contain the direct Pwave. The results are shown
+in Figure \ref{fig_L23D} (left). The L2error reduces with increasing number of time steps $NT=T/\Delta t$ because of the decreasing time step interval $\Delta t$. (The wave propagation time $T=0.07$\,s is held constant).
+At higher orders of the time accuracy ($M$) the error reduces more rapidly proportional to $\Delta t^M$.
+%The observed overall behavior of the convergence of accuracy in the 3D elastic case (Figure \ref{fig_L23D}, left) is quite similar to the convergence obtained for the 1D scheme shown in Figure \ref{fig_L2}.
+%This indicates that the numerical dispersion and dissipation properties derived in chapters
+%\ref{chap1Ddisp} and \ref{chap1Ddiss} for the 1D case are also applicable to Pwaves in 3D media.
+
+In Figure \ref{fig_L23D} (left) we see see that we can reduce the number of time steps that are required to achieve a certain level of accuracy by increasing the temporal accuracy order $M$. For a given error level of $E=0.1\%$ the required number of time steps reduce to approximately $43\%$ and $29\%$ when increasing the temporal order from $M=2$ to $M=3$ and $M=4$, respectively. The corresponding run times decrease less significant because the number of floating point operations also increase with $M$. The observed relation between the number of time steps and the total run time is plotted in Figure \ref{fig_L23D} (right). The computational requirements are also summarized in Table \ref{tab:3D_times}. We see that the corresponding run times reduce only to $58\%$ and $42\%$ for $M=3$ and $M=4$, respectively, which is still a substantial improvement. The downside of these significant run time savings is the boost of the memory requirements which increase to $174\%$ and $220\%$ (Table \ref{tab:3D_times}). Interestingly, the factor for the run time reduction and the increase factor for memory are quite similar for the same accuracy order $M$.
+
+\begin{table}[hb]
+ \begin{center}
+ \begin{tabular}{lllll}
+ & $M=2$ & $M=3$ & $M=4$ & \\\hline
+ Run time & 4349\,s & 2505\,s & 1833\,s & \\
+ & \textbf{100\%} &\textbf{58\%} & \textbf{42\%} & \\\hline
+ Time steps & 8133 & 3505 & 2350 & \\
+ & \textbf{100\%} &\textbf{43\%} & \textbf{29\%} & \\\hline
+ Total memory & 19.9\,GB & 34.7\,GB& 43.8\,GB & \\
+ & \textbf{100\%} & \textbf{174\%} &\textbf{220\%} & \\\hline
+ Error & 0.10013\% &0.10014\% & 0.10013\% & \\
+ \end{tabular}
+ \end{center}
+ \caption{\label{tab:3D_times} Comparison of computational demands for 3D elastic FDTD simulations with the ABS time integrator of accuracy $M$ to achieve the same level of accuracy. The required run time and number of time steps
+ decrease with $M$. At the same time the total requirements of memory increase by approximately the same factor due to the required storage of previously calculated wavefields. }
+\end{table}
+\begin{figure}[ph]
+\centerline{\includegraphics[width=\textwidth]{eps/error_3D_courant}}
+\caption{Left: Relative error of 3D elastic simulations.
+The normalized L2 norm is plotted over the required number of time steps ($NT$) and the courant number. The temporal accuracy orders are $M=2,3,4$. The spatial accuracy is fixed at order $N=8$.
+The number of time steps to achieve a given accuracy of $E=0.1\%$ (horizontal line) reduce to $43\%$ and $29\%$ for the temporal orders $M=3$ and $M=4$, respectively.
+For large $NT$ (small time step intervals) the error converges to the fixed error of the spatial discretization. Right: The required run time of the program as a function of the number of
+time steps. The total run time to achieve an accuracy level of $E=0.1\%$ reduces to $58\%$ and $42\%$ for the temporal orders of $M=3$ and $M=4$, respectively.}
+\label{fig_L23D}
+\end{figure}
\newpage
diff git a/doc/guide_sofi3D/tex/thesis.bib b/doc/guide_sofi3D/tex/thesis.bib
index 673e54da933883f8e60500f8e2edd62ea443fe5d..bcdbdff5be2824dad146b836d3bb7f4188f26166 100644
 a/doc/guide_sofi3D/tex/thesis.bib
+++ b/doc/guide_sofi3D/tex/thesis.bib
@@ 277,11 +277,19 @@ http://www.fzjuelich.de/ias/jsc/EN/Expertise/Supercomputers/JUROPA/JUROPA_node.
YEAR = {1998} }
@inproceedings{bohlen2015higher,
 title={Higher Order FDTD Seismic Modelling Using the Staggered AdamsBashforth Time Integrator},
 author={Bohlen, T and Wittkamp, F},
 booktitle={77th EAGE Conference and Exhibition 2015},
 year={2015}
}
+ Author = {Bohlen, T and Wittkamp, F},
+ Booktitle = {77th EAGE Conference and Exhibition 2015},
+ Title = {Higher Order FDTD Seismic Modelling Using the Staggered AdamsBashforth Time Integrator},
+ Year = {2015}}
+
+@article{bohlen2015higher2,
+ Author = {Bohlen, T and Wittkamp, F},
+ DateAdded = {20150930 10:53:59 +0000},
+ DateModified = {20150930 10:58:00 +0000},
+ Journal = {Geophysical Journal International},
+ Title = {3D viscoelastic FDTD seismic modelling using the staggered AdamsBashforth time integrator},
+ Year = {(submitted)}}
+
@article{bohlen:02,
AUTHOR = {Bohlen, T.},
@@ 546,7 +554,6 @@ http://www.fzjuelich.de/ias/jsc/EN/Expertise/Supercomputers/JUROPA/JUROPA_node.
year={2000},
publisher={SIAM}
}

@article{clayton:80,
AUTHOR = { Clayton, R.W. and Engquist, B.},
TITLE = {Absorbing boundary conditions for waveequation migration},