Skip to content

GitLab

Commit 0e199d69 authored by Florian Wittkamp's avatar Florian Wittkamp
Browse files

Updated Documentation: Adams-Bashforth 

Added a detailed description of the Adams-Bashforth method.
parent ba6e2772
Hide whitespace changes
Inline Side-by-side
AUTHORS
View file @ 0e199d69
......@@ -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
doc/guide_sofi3D/tex/guide_sofi3D.tex
View file @ 0e199d69
......@@ -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 FD-operators}
To achieve a better performance and higher accuracy of the FD simulation higher order temporal FD-operators can be chosen. The order of the temporal integration is increased by the Adams-Bashforth 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 second-order leapfrog method.
Numerical dissipation is introduced by the AB method which is significantly smaller for the AB method of fourth-order accuracy then of third-order accuracy. Moreover the stability limit will be lower for the higher order temporal orders (see \ref{courandt}).\\ \\
A detailed analysis of the Adams-Bashforth method applied to seismic wave simulation with FDTD is given by \cite{bohlen2015higher}.
\subsection{Adams-Bashforth higher-order accurate time integrators}
To improve the precision of the temporal discretization the Adams-Bashforth third-order and fourth-order accurate time integrator can be used.
The staggered Adams-Bashforth method (ABS) is a multi-step 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 second-order leapfrog method. Numerical dissipation is introduced by the ABS method which is significantly smaller for the ABS method of fourth-order accuracy. In 1-D and 3-D simulation experiments the convincing improvements of simulation accuracy of the fourth-order ABS method is verified. In a realistic elastic 3-D 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 Adams-Bashforth method (ABS-method) using the 1-D acoustic wave equation in velocity-stress 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 P-wave 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|^{n-1/2}}{\Delta t}&=& - \pi \frac{\partial v}{\partial x}\bigg |^n + \mathcal{O}(\Delta t^2) \nonumber \\
\frac{v|^n-v|^{n-1}}{\Delta t}&=& -\rho^{-1} \frac{\partial p}{\partial x}\bigg |^{(n-1/2)} + \mathcal{O}(\Delta t^2)
\label{2ndorder}
\end{eqnarray}
This results in the conventional explicit second order accurate time-stepping (leapfrog) scheme
\begin{eqnarray}
p|^{n+1/2}&=&p|^{n-1/2} - \Delta t \pi \frac{\partial v}{\partial x}\bigg |^n + \mathcal{O}(\Delta t^2) \nonumber \\
v|^n&=&v|^{n-1} -\Delta t \rho^{-1} \frac{\partial p}{\partial x}\bigg |^{(n-1/2)} + \mathcal{O}(\Delta t^2)
\label{2ndorderexplicit}
\end{eqnarray}
which requires no additional storage of the wavefield variables $p$ and $v$.
With the ABS-method 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 ABS-method thus requires the storage of previous time levels of spatial derivatives of the pressure $p$ and particle velocity $v$. The time update for the ABS-method thus reads
\begin{eqnarray}
p|^{n+1/2}&=&p|^{n-1/2} - \Delta t \pi \sum_{k=0}^{M-1} a_k \frac{\partial v}{\partial x}\bigg |^{n-k} + \mathcal{O}(\Delta t^M) \nonumber \\
v|^n&=&v|^{n-1} -\Delta t \rho^{-1} \sum_{k=0}^{M-1} a_k \frac{\partial p}{\partial x}\bigg |^{(n-1/2-k)} + \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}{|c|c|c|c|c|}\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 Adams-Bashforth method \protect\cite{ghrist2000staggered}.}
\end{table}
\newpage
\subsubsection{Accuracy in 1-D}
We first compare seismograms calculated with the 1-D ABS-method using equations \ref{absmethod} for different orders of accuracy in time ($M$). We use a 1-D 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 source-receiver 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 (time-shifted 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}, top-left). 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 ABS-method which requires to store $M-1$ 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 L2-misfit 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 L2-misfit 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 1-D 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 multi-step 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 1-D 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 3-D}
In Figure \ref{fig_1dfdseis} we see that the simulation accuracy in 1-D increases with decreasing Courant number and increasing accuracy order. In order to quantify the corresponding numerical simulation error of the 3-D elastic update scheme we calculated the normalized L2-misfit 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 S-velocities, 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 P-wave in a distance of $180$\,m corresponding to 30 dominant wavelength. The size of the model grid is 800x400x400 grid points. The direct P-wave 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 small-scale shared memory cluster.
As a measure of accuracy we use the normalized L2-norm between the numerical and analytical seismograms which only contain the direct P-wave. The results are shown
in Figure \ref{fig_L23D} (left). The L2-error 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 3-D elastic case (Figure \ref{fig_L23D}, left) is quite similar to the convergence obtained for the 1-D 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 1-D case are also applicable to P-waves in 3-D 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}{l|l|l|ll}
& $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 3-D 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 3-D 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
......
doc/guide_sofi3D/tex/thesis.bib
View file @ 0e199d69
......@@ -277,11 +277,19 @@ http://www.fz-juelich.de/ias/jsc/EN/Expertise/Supercomputers/JUROPA/JUROPA_node.
YEAR = {1998} }
@inproceedings{bohlen2015higher,
title={Higher Order FDTD Seismic Modelling Using the Staggered Adams-Bashforth 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 Adams-Bashforth Time Integrator},
Year = {2015}}
@article{bohlen2015higher2,
Author = {Bohlen, T and Wittkamp, F},
Date-Added = {2015-09-30 10:53:59 +0000},
Date-Modified = {2015-09-30 10:58:00 +0000},
Journal = {Geophysical Journal International},
Title = {3-D viscoelastic FDTD seismic modelling using the staggered Adams-Bashforth time integrator},
Year = {(submitted)}}
@article{bohlen:02,
AUTHOR = {Bohlen, T.},
......@@ -546,7 +554,6 @@ http://www.fz-juelich.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 wave-equation migration},
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment