FINITE-DIFFERENCE ALGORITHMS

Explicit finite-difference algorithms map well onto massively parallel computers. At each time step, or depth step, the values of the wavefield is computed by a linear combination of the values of the wavefield at the previous step in the neighboring points of the data grid. Furthermore, the regular data grid used in finite-difference algorithms can be mapped into an hypercube in such a way that the property of locality are preserved. Therefore, only nearest-neighbor communications are required between the processors of an hypercube for executing an explicit finite-difference scheme.

As an example of explicit finite-difference I present a simple subroutine that implements an explicit finite-difference algorithm for modeling acoustic waves. The code shows the subroutine for stepping in time of a simple 2-D second-oder in time and space finite-difference scheme. The values of the wavefield in in the neighboring points are accessed by shifting the data array by use of the Fortran 90 intrinsic function CSHIFT. The most of the neighboring values are stored in the memory of the same processors, and thus no communications is required for accessing them with a CSHIFT. The compiler takes care of issuing the instructions for the interprocessor communications necessary to access the few values that are not local.

	SUBROUTINE TIMESTEP (NT,NX,NY,NSOURCE,IXS,IYS,SOURCEV,
     1	WAVEITM1,WAVEIT,WAVEITP1,C1,C2,C3,C4,C5,C6)
	INTEGER NT,NX,NY,NSOURCE,IXS,IYS

REAL SOURCEV(NSOURCE)
	REAL WAVEITM1(NX,NY),WAVEIT(NX,NY),WAVEITP1(NX,NY)
	REAL C1(NX,NY),C2(NX,NY),C3(NX,NY)
	REAL C4(NX,NY),C5(NX,NY),C6(NX,NY)

C
C time step loop
C

DO IT=1,NT
C
C add point source at location (IXS,IYS)
C

WAVEIT(IXS,IYS)=WAVEIT(IXS,IYS)+SOURCEV(IT)

C
C compute finite-difference star with spatially varying coefficients
C C1,C2,C3,C4,C5
C

WAVEITP1 = C1*CSHIFT(WAVEIT,DIM=1,SHIFT=-1)
     1			+ C2*WAVEIT
     2     		+ C3*CSHIFT(WAVEIT,DIM=1,SHIFT= 1)
     3    		+ C4*CSHIFT(WAVEIT,DIM=2,SHIFT=-1)
     4   		+ C5*CSHIFT(WAVEIT,DIM=2,SHIFT= 1)

C
C compute next time step
C
		WAVEITP1=WAVEITP1+C6*WAVEITM1

C
C update value for previous time step
C

WAVEITM1=WAVEIT
		WAVEIT=WAVEITP1

END DO

RETURN
	END

The efficient treatment of boundary conditions can be a problem on SIMD machines, because the computations performed in the boundary regions may be different than the computations executed in the interior region. In this case, the computations in the boundary must be performed after the computations in the interior, and not in parallel . Often this problem can be solved by applying the same finite-difference star in the boundary region as in the interior region, but with different coefficients. Changing the coefficients makes the finite-difference equation solved on the boundary different than the equation solved in the interior, without requiring a different differencing star. In the example shown before the variability of the coefficients have been exploited for implementing free-surface boundary conditions at the top of the computational domain, and absorbing boundary conditions on the other three sides. Figure  shows a snapshot of the wavefield generated using this subroutine.

The efficient implementation of implicit finite-difference algorithms is not always as straightforward as explicit schemes, because implicit algorithms require the solution of a sparse system of linear equations. The linear systems can be solved by iterative algorithm like conjugate gradient (Nichols, 1991) The iterative solution are usually easy to parallelize because they require computations similar to the one required by explicit finite differences. On the contrary, the commonly used algorithms for direct solution of sparse linear systems are recursive. Important example in Geophysics is the solution of the tridiagonal system in finite-difference migration, that creates problems for vector computers as well as parallel computers. On a parallel computer the tridiagonal systems can be solved with methods similar to the ones used on vector computers, like cyclic reduction (Johnsson, 1987, Cole 1991). Or in some cases, as for example for finite-difference prestack migration, the particular structure of the problem can be exploited for devising a non-recursive parallel algorithm (Moorhead and Biondi, 1991).

 

fdfig
Figure 2 Snapshot of the wavefield captured from the CM FrameBuffer. Notice the absorbing boundary conditions on the sides and the free surface at the top.