[ < ] | [ > ] | [ << ] | [ Up ] | [ >> ] | [Top] | [Contents] | [Index] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
central12(A,dimension)
1st derivative, 2nd order accurate. Factor: 2h
-1 | 0 | 1 | |
-1 | 1 |
central22(A,dimension)
2nd derivative, 2nd order accurate. Factor: h^2
-1 | 0 | 1 | |
1 | -2 | 1 |
central32(A,dimension)
3rd derivative, 2nd order accurate. Factor: 2h^3
-2 | -1 | 0 | 1 | 2 | |
-1 | 2 | -2 | 1 |
central42(A,dimension)
4th derivative, 2nd order accurate. Factor: h^4
-2 | -1 | 0 | 1 | 2 | |
1 | -4 | 6 | -4 | 1 |
central14(A,dimension)
1st derivative, 4th order accurate. Factor: 12h
-2 | -1 | 0 | 1 | 2 | |
1 | -8 | 8 | -1 |
central24(A,dimension)
2nd derivative, 4th order accurate. Factor: 12h^2
-2 | -1 | 0 | 1 | 2 | |
-1 | 16 | -30 | 16 | -1 |
central34(A,dimension)
3rd derivative, 4th order accurate. Factor: 8h^3
-2 | -1 | 0 | 1 | 2 | |
-8 | 13 | -13 | 8 |
central44(A,dimension)
4th derivative, 4th order accurate. Factor: 6h^4
-2 | -1 | 0 | 1 | 2 | |
12 | -39 | 56 | -39 | 12 |
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] | [ ? ] |
forward11(A,dimension)
1st derivative, 1st order accurate. Factor: h
0 | 1 | |
-1 | 1 |
forward21(A,dimension)
2nd derivative, 1st order accurate. Factor: h^2
0 | 1 | 2 | |
1 | -2 | 1 |
forward31(A,dimension)
3rd derivative, 1st order accurate. Factor: h^3
0 | 1 | 2 | 3 | |
-1 | 3 | -3 | 1 |
forward41(A,dimension)
4th derivative, 1st order accurate. Factor: h^4
0 | 1 | 2 | 3 | 4 | |
1 | -4 | 6 | -4 | 1 |
forward12(A,dimension)
1st derivative, 2nd order accurate. Factor: 2h
0 | 1 | 2 | |
-3 | 4 | -1 |
forward22(A,dimension)
2nd derivative, 2nd order accurate. Factor: h^2
0 | 1 | 2 | 3 | |
2 | -5 | 4 | -1 |
forward32(A,dimension)
3rd derivative, 2nd order accurate. Factor: 2h^3
0 | 1 | 2 | 3 | 4 | |
-5 | 18 | -24 | 14 | -3 |
forward42(A,dimension)
4th derivative, 2nd order accurate. Factor: h^4
0 | 1 | 2 | 3 | 4 | 5 | |
3 | -14 | 26 | -24 | 11 | -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] | [ ? ] |
backward11(A,dimension)
1st derivative, 1st order accurate. Factor: h
-1 | 0 | |
-1 | 1 |
backward21(A,dimension)
2nd derivative, 1st order accurate. Factor: h^2
-2 | -1 | 0 | |
1 | -2 | 1 |
backward31(A,dimension)
3rd derivative, 1st order accurate. Factor: h^3
-3 | -2 | -1 | 0 | |
-1 | 3 | -3 | 1 |
backward41(A,dimension)
4th derivative, 1st order accurate. Factor: h^4
-4 | -3 | -2 | -1 | 0 | |
1 | -4 | 6 | -4 | 1 |
backward12(A,dimension)
1st derivative, 2nd order accurate. Factor: 2h
-2 | -1 | 0 | |
1 | -4 | 3 |
backward22(A,dimension)
2nd derivative, 2nd order accurate. Factor: h^2
-3 | -2 | -1 | 0 | |
-1 | 4 | -5 | 2 |
backward32(A,dimension)
3rd derivative, 2nd order accurate. Factor: 2h^3
-4 | -3 | -2 | -1 | 0 | |
3 | -14 | 24 | -18 | 5 |
backward42(A,dimension)
4th derivative, 2nd order accurate. Factor: h^4
-5 | -4 | -3 | -2 | -1 | 0 | |
-2 | 11 | -24 | 26 | -14 | 3 |
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] | [ ? ] |
Laplacian2D(A)
2nd order accurate, 2-dimensional laplacian. Factor: h^2
-1 | 0 | 1 | |
-1 | 1 | ||
0 | 1 | -4 | 1 |
1 | 1 |
Laplacian3D(A)
2nd order accurate, 3-dimensional laplacian. Factor: h^2
Laplacian2D4(A)
4th order accurate, 2-dimensional laplacian. Factor: 12h^2
-2 | -1 | 0 | 1 | 2 | |
-2 | -1 | ||||
-1 | 16 | ||||
0 | -1 | 16 | -60 | 16 | -1 |
1 | 16 | ||||
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] | [ ? ] |
These return TinyVector
s 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] | [ ? ] |
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] | [ ? ] |
There are also grad-squared operators, which return TinyVector
s 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] | [ ? ] |
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 TinyVector
s 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] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
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] | [ ? ] |
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.