RTM using effective boundary saving: A staggered grid GPU implementation

Effective boundary for regular grid finite difference

Assume $2N$ -th order finite difference scheme is applied. The Laplacian operator is specified by

$$\begin{array}{rl} \nabla^2 p^{k}&=\partial_{xx}p^{k}+\partial... ...frac{1}{\Delta x^2}\sum_{i=-N}^Nc_i p^k[ix+i][iz]\\ \end{array}$$ (6)

where $c_i$ is given by Table 1, see a detailed derivation in Fornberg (1988). The Laplacian operator has $x$ and $z$ with same finite difference structure. For $x$ dimension only, the second derivative of order $2N$ requires at least $N$ points in the boundary zone, as illustrated by Figure 2. In 2-D case, the required boundary zone has been plotted in Figure 3a. Note that four corners in $B_1B_2B_3B_4$ in Figure 1 are not needed. This is exactly the boundary saving scheme proposed by Dussaud et al. (2008).

Table 1: Finite difference coefficients for regular grid (Order-$2N$ )

$i$	-4	-3	-2	-1	0	1	2	3	4
$N=1$				1	-2	1
$N=2$			-1/12	4/3	-5/2	4/3	-1/12
$N=3$		1/90	-3/20	3/2	-49/18	3/2	-3/20	1/90
$N=4$	-1/560	8/315	-1/5	8/5	-205/72	8/5	-1/5	8/315	-1/560

Keep in mind that we only need to guarantee the correctness of the wavefield in the original model zone $A_1A_2A_3A_4$ . However, the saved wavefield in $A_1A_2A_3A_4\backslash B_1B_2B_3B_4$ is also correct. Is it possible to further shrink it to reduce number of points for saving? The answer is true. Our solution is: saving the inner $N$ layers on each side neighboring the boundary $A_1A_2A_3A_4\backslash D_1D_2D_3D_4$ , as shown in Figure 3b. We call it the effective boundary for regular finite difference scheme.

After $nt$ steps of forward modeling, we begin our backward propagation with the last 2 wavefield snap $p^{nt}$ and $p^{nt-1}$ and saved effective boundaries in $A_1A_2A_3A_4\backslash D_1D_2D_3D_4$ . At that moment, the wavefield is correct for every grid point. (Of course, the correctness of the wavefield in $A_1A_2A_3A_4$ is guaranteed.) At time $k$ , we assume the wavefield in $A_1A_2A_3A_4$ is correct. One step of backward propagation means $A_1A_2A_3A_4$ is shrunk to $D_1D_2D_3D_4$ . In other words, the wavefield in $D_1D_2D_3D_4$ is correctly reconstructed. Then we load the saved effective boundary of time $k$ to overwrite the area $A_1A_2A_3A_4\backslash D_1D_2D_3D_4$ . Again, all points of the wavefield in $A_1A_2A_3A_4$ are correct. We repeat this overwriting and computing process from one time step to another ( $k\rightarrow k-1$ ), in reverse time order. The wavefield in the boundary $C_1C_2C_3C_4\backslash A_1A_2A_3A_4$ may be incorrect because the points here are neither saved nor correctly reconstructed from the previous step.

fig2
Figure 2. 1-D schematic plot of required points in regular grid for boundary saving. Computing the laplacian needs $N$ points in the extended boundary zone, the rest $N+1$ points in the inner model grid. $N$ points is required for boundary saving.

fig3
Figure 3. A 2-D sketch of required points for boundary saving for regular grid finite difference: (a) The scheme proposed by Dussaud et al. (2008) (red zone). (b) Proposed effective boundary saving scheme (gray zone).

RTM using effective boundary saving: A staggered grid GPU implementation