arrays-stencils.texi

A C++ class library for scientific computing

@cindex stencil objects
@cindex Array stencils

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

@section 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:

@example
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);
@end example

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

@example
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);
@end example

@section Declaring stencil objects
@cindex stencil objects declaring

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 @code{BZ_DECLARE_STENCIL1},
@code{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 @code{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 @code{BZ_END_STENCIL} must appear
(or the macro @code{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:

@findex BZ_DECLARE_STENCIL

@example
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
@end example

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

You can use the notation @code{A(i,j,k)} to read the element at an offset
@code{(i,j,k)} from the current element.  If you omit the parentheses
(i.e.@: as in ``@code{A}'' then the current element is read.

You can invoke @emph{stencil operators} which calculate finite differences
and laplacians.

@section Automatic determination of stencil extent

In stencil declarations such as

@example
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
@end example

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
@code{BZ_END_STENCIL_WITH_SHAPE}:

@example
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))
@end example

The parameters of this macro are: a @code{TinyVector} (constructed by the
@code{shape()} function) containing the lower bounds of the stencil offsets,
and a @code{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:

@example
      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
@end example

@section Stencil operators
@cindex 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 @code{h}.  A 2nd-order accurate operator has error term 
@math{O(h^2)}; a 4th-order accurate operator has error term @math{O(h^4)}.

All of the stencils have factors associated with them.  For example, the
@code{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 @code{central12(A,firstDim)}@math{/(2h)}.
Typically when designing stencils, one factors out all of the @math{h} terms
for efficiency.

The factor terms always consist of an integer multiplier (often 1) and a
power of @math{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 @code{n} appended to the name, for
example @code{central12n} is the normalized version of @code{central12}, and
has factor @math{h} instead of @math{2h}.

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

@subsection Central differences
@cindex central differences

@table @code
@item central12(A,dimension)
1st derivative, 2nd order accurate.  Factor: @math{2h}
@include stencils/central12.texi

@item central22(A,dimension)
2nd derivative, 2nd order accurate.  Factor: @math{h^2}
@include stencils/central22.texi

@item central32(A,dimension)
3rd derivative, 2nd order accurate.  Factor: @math{2h^3}
@include stencils/central32.texi

@item central42(A,dimension)
4th derivative, 2nd order accurate.  Factor: @math{h^4}
@include stencils/central42.texi

@item central14(A,dimension)
1st derivative, 4th order accurate.  Factor: @math{12h}
@include stencils/central14.texi

@item central24(A,dimension)
2nd derivative, 4th order accurate.  Factor: @math{12h^2}
@include stencils/central24.texi

@item central34(A,dimension)
3rd derivative, 4th order accurate.  Factor: @math{8h^3}
@include stencils/central34.texi

@item central44(A,dimension)
4th derivative, 4th order accurate.  Factor: @math{6h^4}
@include stencils/central44.texi
@end table

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

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

@subsection Forward differences
@cindex forward differences

@table @code
@item forward11(A,dimension)
1st derivative, 1st order accurate.  Factor: @math{h}
@include stencils/forward11.texi

@item forward21(A,dimension)
2nd derivative, 1st order accurate.  Factor: @math{h^2}
@include stencils/forward21.texi

@item forward31(A,dimension)
3rd derivative, 1st order accurate.  Factor: @math{h^3}
@include stencils/forward31.texi

@item forward41(A,dimension)
4th derivative, 1st order accurate.  Factor: @math{h^4}
@include stencils/forward41.texi

@item forward12(A,dimension)
1st derivative, 2nd order accurate.  Factor: @math{2h}
@include stencils/forward12.texi

@item forward22(A,dimension)
2nd derivative, 2nd order accurate.  Factor: @math{h^2}
@include stencils/forward22.texi

@item forward32(A,dimension)
3rd derivative, 2nd order accurate.  Factor: @math{2h^3}
@include stencils/forward32.texi

@item forward42(A,dimension)
4th derivative, 2nd order accurate.  Factor: @math{h^4}
@include stencils/forward42.texi
@end table

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

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

@subsection Backward differences
@cindex backward differences

@table @code
@item backward11(A,dimension)
1st derivative, 1st order accurate.  Factor: @math{h}
@include stencils/backward11.texi

@item backward21(A,dimension)
2nd derivative, 1st order accurate.  Factor: @math{h^2}
@include stencils/backward21.texi

@item backward31(A,dimension)
3rd derivative, 1st order accurate.  Factor: @math{h^3}
@include stencils/backward31.texi

@item backward41(A,dimension)
4th derivative, 1st order accurate.  Factor: @math{h^4}
@include stencils/backward41.texi

@item backward12(A,dimension)
1st derivative, 2nd order accurate.  Factor: @math{2h}
@include stencils/backward12.texi

@item backward22(A,dimension)
2nd derivative, 2nd order accurate.  Factor: @math{h^2}
@include stencils/backward22.texi

@item backward32(A,dimension)
3rd derivative, 2nd order accurate.  Factor: @math{2h^3}
@include stencils/backward32.texi