THE ADAPTIVE STEP IMPLEMENTATION

The finite difference algorithm described by Van Trier and Symes (1991) is based on the eikonal equation in cartesian coordinates

 

u²+v²=s² (1)

where

$$\left \{ \begin{array} {l} u={\partial{t(x,z)} \over \partia... ...\\ \\ v={\partial{t(x,z)} \over \partial z} \end{array} \right.$$

and where s(x,z) is the 2-dimensional slowness model and t(x,z) is the travel time field. In cylindrical coordinates the eikonal equation becomes  

$${u^2 \over r^2}+v^2=s^2$$ (2)

where

$$\left \{ \begin{array} {l} u={\partial{t(r,\theta)} \over \p... ...{\partial{t(r,\theta)} \over \partial r} \\ \end{array} \right.$$

The traveltime is propagated one depth step at a time in cartesian coordinates, or one radial step at a time in cylindrical coordinates. The values of $u={\Delta{t} \over \Delta x}$ and $v={{\Delta{t}} \over {\Delta z}}$ are calculated on each computational front. The domain of dependence is given by the length ratio of the two vectors $\vec{u}$ and $\vec{v}$, as seen in Figure 1. The fixed grid algorithm is faster than the adaptive algorithm, sometimes by a factor of 10-20, depending on the velocity model considered. However numerical instabilities like the one in Figure 2 can not be predicted. The slower adaptive step algorithm will offer robust results.

 

TravelTimes
Figure 2 Two cases for a constant gradient velocity model. a) Numerical instabilities appear in the fixed grid algorithm. b) Given the same initial grid the adaptive step algorithm provides the correct solution.

My initial adaptive step calculation in polar coordinates, based on the geometrical arguments presented in Figure 1, is listed below.

The calculations are done for each new radial computational front. The variable ddr is the initial radial increment which is modified adaptively. The arrays v and u store the values of the functions defined in equation (2). The minimum angle for the domain of dependence is calculated for each computational front and the minimum radial increment is used to advance the front.

A more elegant approach was suggested by Dave Hale. It monitors the same domain of dependence condition avoiding all the time consuming trigonometric relations. It also avoids the division by zero that can occur in either the case of u=0 or v=0. Essentially, the trick is to check for the ratio between v and $u \over r$ in the form v * r <= TINY * u which increases the dynamic range of the variables compared.