4.2 Declaring stencil objects

A stencil declaration may not be inside a function. It can appear inside a class declaration (in which case the stencil object is a nested type).

Stencil objects are declared using the macros BZ_DECLARE_STENCIL1, BZ_DECLARE_STENCIL2, etc. The number suffix is how many arrays are involved in the stencil (in the above example, 4 arrays– P1, P2, P3, c – are used, so the macro BZ_DECLARE_STENCIL4 is invoked).

The first argument is a name for the stencil object. Subsequent arguments are names for the arrays on which the stencil operates.

After the stencil declaration, the macro BZ_END_STENCIL must appear (or the macro BZ_END_STENCIL_WITH_SHAPE, described in the next section).

In between the two macros, you can have multiple assignment statements, if/else/elseif constructs, function calls, loops, etc.

Here are some simple examples:

BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL

BZ_DECLARE_STENCIL4(acoustic2D,P1,P2,P3,c)
  A = 2 * P2 + c * (-4 * P2(0,0) + P2(0,1) + P2(0,-1) + P2(1,0) + P2(-1,0))
      - P1;
BZ_END_STENCIL

BZ_DECLARE_STENCIL8(prop2D,E1,E2,E3,M1,M2,M3,cE,cM)
  E3 = 2 * E2 + cE * Laplacian2D(E2) - E1;
  M3 = 2 * M2 + cM * Laplacian2D(M2) - M1;
BZ_END_STENCIL

BZ_DECLARE_STENCIL3(smooth2Db,A,B,c)
  if ((c > 0.0) && (c < 1.0))
    A = c * (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0
      + (1-c)*B;
  else
    A = 0;
BZ_END_STENCIL

Currently, a stencil can take up to 11 array parameters.

You can use the notation A(i,j,k) to read the element at an offset (i,j,k) from the current element. If you omit the parentheses (i.e. as in “A” then the current element is read.

You can invoke stencil operators which calculate finite differences and laplacians.

4.3 Automatic determination of stencil extent

In stencil declarations such as

BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL

Blitz++ will try to automatically determine the spatial extent of the stencil. This will usually work for stencils defined on integer or float arrays. However, the mechanism does not work well for complex-valued arrays, or arrays of user-defined types. If you get a peculiar error when you try to use a stencil, you probably need to tell Blitz++ the special extent of the stencil manually.

You do this by ending a stencil declaration with BZ_END_STENCIL_WITH_SHAPE:

BZ_DECLARE_STENCIL2(smooth2D,A,B)
  A = (B(0,0) + B(0,1) + B(0,-1) + B(1,0) + B(-1,0)) / 5.0;
BZ_END_STENCIL_WITH_SHAPE(shape(-1,-1),shape(+1,+1))

The parameters of this macro are: a TinyVector (constructed by the shape() function) containing the lower bounds of the stencil offsets, and a TinyVector containing the upper bounds. You can determine this by looking at the the terms in the stencil and finding the minimum and maximum value of each index:

      A = (B(0,  0) 
         + B(0, +1)
         + B(0, -1)
         + B(+1, 0)
         + B(-1, 0)) / 5.0;
           --------
min indices  -1, -1
max indices  +1, +1