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

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,IYSREAL 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).

Figure 2

12/18/1997