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