Blitz logo

Blitz Support :

From: Faheem Mitha (faheem_at_[hidden])
Date: 2005-05-30 16:05:53


Hi,

In the manual, consider the following section. I don't understand how the
former expression

P3(I,J,K) = (2-6*c(I,J,K)) * P2(I,J,K)
             + c(I,J,K)*(P2(I-1,J,K) + P2(I+1,J,K) + P2(I,J-1,K) + P2(I,J+1,K)
             + P2(I,J,K-1) + P2(I,J,K+1)) - P1(I,J,K);

is equivalent to the latter

applyStencil(acoustic3D_stencil(), P1, P2, P3, c);

Should this perhaps be

applyStencil(acoustic3D_stencil(), P1(I,J,K), P2(I,J,K), P3(I,J,K), c);

instead? Otherwise, how does the stencil know how to restrict to the
appropriate range?

To put it more crudely, how can the two expressions be equivalent,
when the former is written in terms of the Range objects I, J, K, and
the latter is not?

Also, I'm unclear on the difference between BZ_DECLARE_STENCIL and
BZ_DECLARE_STENCIL_OPERATOR. The latter is only defined for up to 3
arguments, but I don't see any difference except in the way they are
called. Perhaps only stencil operators can be called inside a stencil
declaration?

Thanks in advance for any clarification.

                                                                  Faheem.

**************************************************************************
4.1 Motivation: a nicer notation for stencils

Suppose we wanted to implement the 3-D acoustic wave equation using finite
differencing. Here is how a single iteration would look using subarray
syntax:

Range I(1,N-2), J(1,N-2), K(1,N-2);

P3(I,J,K) = (2-6*c(I,J,K)) * P2(I,J,K)
             + c(I,J,K)*(P2(I-1,J,K) + P2(I+1,J,K) + P2(I,J-1,K) + P2(I,J+1,K)
             + P2(I,J,K-1) + P2(I,J,K+1)) - P1(I,J,K);

This syntax is a bit klunky. With stencil objects, the implementation
becomes:

BZ_DECLARE_STENCIL4(acoustic3D_stencil,P1,P2,P3,c)
   P3 = 2 * P2 + c * Laplacian3D(P2) - P1;
BZ_END_STENCIL
   .
   .
applyStencil(acoustic3D_stencil(), P1, P2, P3, c);