[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4. Stencils

Blitz++ provides an implementation of stencil objects which is currently experimental. This means that the exact details of how they are declared and used may change in future releases. Use at your own risk.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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);

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

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

[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4 Stencil operators

This section lists all the stencil operators provided by Blitz++. They assume that an array represents evenly spaced data points separated by a distance of h. A 2nd-order accurate operator has error term O(h^2); a 4th-order accurate operator has error term O(h^4).

All of the stencils have factors associated with them. For example, the central12 operator is a discrete first derivative which is 2nd-order accurate. Its factor is 2h; this means that to get the first derivative of an array A, you need to use central12(A,firstDim)/(2h). Typically when designing stencils, one factors out all of the h terms for efficiency.

The factor terms always consist of an integer multiplier (often 1) and a power of h. For ease of use, all of the operators listed below are provided in a second "normalized" version in which the integer multiplier is 1. The normalized versions have an n appended to the name, for example central12n is the normalized version of central12, and has factor h instead of 2h.

These operators are defined in blitz/array/stencilops.h if you wish to see the implementation.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.1 Central differences

central12(A,dimension)

1st derivative, 2nd order accurate. Factor: 2h

-101
-11
central22(A,dimension)

2nd derivative, 2nd order accurate. Factor: h^2

-101
1-21
central32(A,dimension)

3rd derivative, 2nd order accurate. Factor: 2h^3

-2-1012
-12-21
central42(A,dimension)

4th derivative, 2nd order accurate. Factor: h^4

-2-1012
1-46-41
central14(A,dimension)

1st derivative, 4th order accurate. Factor: 12h

-2-1012
1-88-1
central24(A,dimension)

2nd derivative, 4th order accurate. Factor: 12h^2

-2-1012
-116-3016-1
central34(A,dimension)

3rd derivative, 4th order accurate. Factor: 8h^3

-2-1012
-813-138
central44(A,dimension)

4th derivative, 4th order accurate. Factor: 6h^4

-2-1012
12-3956-3912

Note that the above are available in normalized versions central12n, central22n, ..., central44n which have factors of h, h^2, h^3, or h^4 as appropriate.

These are available in multicomponent versions: for example, central12(A,component,dimension) gives the central12 operator for the specified component (Components are numbered 0, 1, ... N-1).


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.2 Forward differences

forward11(A,dimension)

1st derivative, 1st order accurate. Factor: h

01
-11
forward21(A,dimension)

2nd derivative, 1st order accurate. Factor: h^2

012
1-21
forward31(A,dimension)

3rd derivative, 1st order accurate. Factor: h^3

0123
-13-31
forward41(A,dimension)

4th derivative, 1st order accurate. Factor: h^4

01234
1-46-41
forward12(A,dimension)

1st derivative, 2nd order accurate. Factor: 2h

012
-34-1
forward22(A,dimension)

2nd derivative, 2nd order accurate. Factor: h^2

0123
2-54-1
forward32(A,dimension)

3rd derivative, 2nd order accurate. Factor: 2h^3

01234
-518-2414-3
forward42(A,dimension)

4th derivative, 2nd order accurate. Factor: h^4

012345
3-1426-2411-2

Note that the above are available in normalized versions forward11n, forward21n, ..., forward42n which have factors of h, h^2, h^3, or h^4 as appropriate.

These are available in multicomponent versions: for example, forward11(A,component,dimension) gives the forward11 operator for the specified component (Components are numbered 0, 1, ... N-1).


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.3 Backward differences

backward11(A,dimension)

1st derivative, 1st order accurate. Factor: h

-10
-11
backward21(A,dimension)

2nd derivative, 1st order accurate. Factor: h^2

-2-10
1-21
backward31(A,dimension)

3rd derivative, 1st order accurate. Factor: h^3

-3-2-10
-13-31
backward41(A,dimension)

4th derivative, 1st order accurate. Factor: h^4

-4-3-2-10
1-46-41
backward12(A,dimension)

1st derivative, 2nd order accurate. Factor: 2h

-2-10
1-43
backward22(A,dimension)

2nd derivative, 2nd order accurate. Factor: h^2

-3-2-10
-14-52
backward32(A,dimension)

3rd derivative, 2nd order accurate. Factor: 2h^3

-4-3-2-10
3-1424-185
backward42(A,dimension)

4th derivative, 2nd order accurate. Factor: h^4

-5-4-3-2-10
-211-2426-143

Note that the above are available in normalized versions backward11n, backward21n, ..., backward42n which have factors of h, h^2, h^3, or h^4 as appropriate.

These are available in multicomponent versions: for example, backward42(A,component,dimension) gives the backward42 operator for the specified component (Components are numbered 0, 1, ... N-1).


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.4 Laplacian (\nabla

Laplacian2D(A)

2nd order accurate, 2-dimensional laplacian. Factor: h^2

-101
-11
01-41
11
Laplacian3D(A)

2nd order accurate, 3-dimensional laplacian. Factor: h^2

Laplacian2D4(A)

4th order accurate, 2-dimensional laplacian. Factor: 12h^2

-2-1012
-2-1
-116
0-116-6016-1
116
2-1
Laplacian3D4(A)

4th order accurate, 3-dimensional laplacian. Factor: 12h^2

Note that the above are available in normalized versions Laplacian2D4n, Laplacian3D4n which have factors h^2.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.5 Gradient (\nabla

These return TinyVectors of the appropriate numeric type and length:

grad2D(A)

2nd order, 2-dimensional gradient (vector of first derivatives), generated using the central12 operator. Factor: 2h

grad2D4(A)

4th order, 2-dimensional gradient, using central14 operator. Factor: 12h

grad3D(A)

2nd order, 3-dimensional gradient, using central12 operator. Factor: 2h

grad3D4(A)

4th order, 3-dimensional gradient, using central14 operator. Factor: 12h

These are available in normalized versions grad2Dn, grad2D4n, grad3Dn and grad3D4n which have factors h.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.6 Jacobian operators

The Jacobian operators are defined over 3D vector fields only (e.g. Array<TinyVector<double,3>,3>). They return a TinyMatrix<T,3,3> where T is the numeric type of the vector field.

Jacobian3D(A)

2nd order, 3-dimensional Jacobian using the central12 operator. Factor: 2h.

Jacobian3D4(A)

4th order, 3-dimensional Jacobian using the central14 operator. Factor: 12h.

These are also available in normalized versions Jacobian3Dn and Jacobain3D4n which have factors h.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.7 Grad-squared operators

There are also grad-squared operators, which return TinyVectors of second derivatives:

gradSqr2D(A)

2nd order, 2-dimensional grad-squared (vector of second derivatives), generated using the central22 operator. Factor: h^2

gradSqr2D4(A)

4th order, 2-dimensional grad-squared, using central24 operator. Factor: 12h^2

gradSqr3D(A)

2nd order, 3-dimensional grad-squared, using the central22 operator. Factor: h^2

gradSqr3D4(A)

4th order, 3-dimensional grad-squared, using central24 operator. Factor: 12h^2

Note that the above are available in normalized versions gradSqr2Dn, gradSqr2D4n, gradSqr3Dn, gradSqr3D4n which have factors h^2.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.8 Curl (\nabla

These curl operators return scalar values:

curl(Vx,Vy)

2nd order curl operator using the central12 operator. Factor: 2h

curl4(Vx,Vy)

4th order curl operator using the central14 operator. Factor: 12h

curl2D(V)

2nd order curl operator on a 2D vector field (e.g. Array<TinyVector<float,2>,2>), using the central12 operator. Factor: 2h

curl2D4(V)

4th order curl operator on a 2D vector field, using the central12 operator. Factor: 12h

Available in normalized forms curln, curl4n, curl2Dn, curl2D4n.

These curl operators return three-dimensional TinyVectors of the appropriate numeric type:

curl(Vx,Vy,Vz)

2nd order curl operator using the central12 operator. Factor: 2h

curl4(Vx,Vy,Vz)

4th order curl operator using the central14 operator. Factor: 12h

curl(V)

2nd order curl operator on a 3D vector field (e.g. Array<TinyVector<double,3>,3>, using the central12 operator. Factor: 2h

curl4(V)

4th order curl operator on a 3D vector field, using the central14 operator. Factor: 12h

Note that the above are available in normalized versions curln and curl4n, which have factors of h.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.9 Divergence (\nabla

The divergence operators return a scalar value.

div(Vx,Vy)

2nd order div operator using the central12 operator. Factor: 2h

div4(Vx,Vy)

4th order div operator using the central14 operator. Factor: 12h

div(Vx,Vy,Vz)

2nd order div operator using the central12 operator. Factor: 2h

div4(Vx,Vy,Vz)

4th order div operator using the central14 operator. Factor: 12h

div2D(V)

2nd order div operator on a 2D vector field, using the central12 operator. Factor: 2h

div2D4(V)

2nd order div operator on a 2D vector field, using the central14 operator. Factor: 12h

div3D(V)

2nd order div operator on a 3D vector field, using the central12 operator. Factor: 2h

div3D4(V)

2nd order div operator on a 3D vector field using the central14 operator. Factor: 12h

These are available in normalized versions divn, div4n, div2Dn, div2D4n, div3Dn, and div3D4n which have factors of h.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.4.10 Mixed partial derivatives

mixed22(A,dim1,dim2)

2nd order accurate, 2nd mixed partial derivative. Factor: 4h^2

mixed24(A,dim1,dim2)

4th order accurate, 2nd mixed partial derivative. Factor: 144h^2

There are also normalized versions of the above, mixed22n and mixed24n which have factors h^2.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.5 Declaring your own stencil operators

You can declare your own stencil operators using the macro BZ_DECLARE_STENCIL_OPERATOR1. For example, here is the declaration of Laplacian2D:

 
BZ_DECLARE_STENCIL_OPERATOR1(Laplacian2D, A)
    return -4*A(0,0) + A(-1,0) + A(1,0) + A(0,-1) + A(0,1);
BZ_END_STENCIL_OPERATOR

To declare a stencil operator on 3 operands, use the macro BZ_DECLARE_STENCIL_OPERATOR3. Here is the declaration of div:

 
BZ_DECLARE_STENCIL_OPERATOR3(div,vx,vy,vz)
  return central12(vx,firstDim) + central12(vy,secondDim)
    + central12(vz,thirdDim);
BZ_END_STENCIL_OPERATOR

The macros aren't magical; they just declare an inline template function with the names and arguments you specify. For example, the declaration of div could also be written

 
template<class T>                              
inline typename T::T_numtype div(T& vx, T& vy, T& vz)   
{
  return central12(vx,firstDim) + central12(vy,secondDim)
                                + central12(vz,thirdDim);
}

The template parameter T is an iterator type for arrays.

You are encouraged to use the macros when possible, because it is possible the implementation could be changed in the future.

To declare a difference operator, use this syntax:

 
BZ_DECLARE_DIFF(central12,A) {
  return A.shift(1,dim) - A.shift(-1,dim);
}

The method shift(offset,dim) retrieves the element at offset in dimension dim.

Stencil operator declarations cannot occur inside a function. If declared inside a class, they are scoped by the class.


[ < ] [ > ]   [ << ] [ Up ] [ >> ]         [Top] [Contents] [Index] [ ? ]

4.6 Applying a stencil

The syntax for applying a stencil is:

 
applyStencil(stencilname(),A,B,C...,F);

Where stencilname is the name of the stencil, and A,B,C,...,F are the arrays on which the stencil operates.

For examples, see `examples/stencil.cpp' and `examples/stencil2.cpp'.

Blitz++ interrogates the stencil object to find out how large its footprint is. It only applies the stencil over the region of the arrays where it won't overrun the boundaries.


[ << ] [ >> ]           [Top] [Contents] [Index] [ ? ]

This document was generated by Charlie & on April, 4 2005 using texi2html 1.76.