|
|
Blitz Devel : |
From: Todd Veldhuizen (tveldhui_at_[hidden])
Date: 1998-07-22 17:00:30
I've implemented a 4th-order computational fluid dynamics solver
using the new multicomponent arrays and stencil objects in blitz.
The source code is available as
http://oonumerics.org/blitz/examples/Blitz++/cfd.cpp
Documentation for stencil objects:
http://oonumerics.org/blitz/manual/blitz02.html#l63
For multicomponent arrays:
http://oonumerics.org/blitz/manual/blitz02.html#l70
Some things work fairly well, others don't. I'm interested
in feedback, particularly if you are a PDE/finite difference
person.
Good things:
- multicomponent arrays make it easy to work with vector fields,
e.g.
typedef TinyVector<double,3> vectorField;
Array<vectorField,3> force;
const int x=0, y=1, z=2;
force[x] = 0;
force[y] = 0;
force[z] = gravity;
- stencil operators take away the drudgery of writing out
big stencils. e.g. the advection calculation looks like this:
advect = product(grad3DVec4(V,geom), *V);
where V is the velocity field, and "geom" is the geometry object.
The stencil operator grad3DVec4 takes a 4th-order accurate gradient
of a vector field, and is really a 36-point stencil or somesuch.
It returns a 3x3 matrix of derivatives, which then multiplies the
velocity field.
- stencil objects let you pass array expressions to functions. So
the conjugate gradient solver takes a stencil object, the arrays,
and the boundary conditions as parameters. This makes it nicely
generic: you can solve any system of equations on any shape/rank
of array, with any BCs.
- stencil objects solve the "Quinlan problem" of getting 1 new pointer
for every stencil point, and resulting register spillage.
But there is also a lot of ugliness:
- arrays, geometry, and boundary conditions are separate objects.
This means that stencil operators need to be given the geometry too
e.g. grad3DVec4(V,geom) instead of just grad3DVec4(V). On the plus
side, geometry objects can be shared among multiple arrays, and
stencil operators can have different implementations for different
geometries.
- declaring a stencil object for each expression is really messy.
Need to "expression templatize" these beasts.
- currently can't have scalar parameters for a stencil object; this
forces you to make all the scalars global variables (super ugh).
There is a workaround, but I haven't implemented it yet.
- blitz isn't completely happy about arrays of vectors; in some
cases you have to stick in "*"'s to explicitly dereference an
iterator, e.g.
advect = product(grad3DVec4(V,geom), *V);
^^^
It's a bit unpredictable where you need the *'s and where you don't.
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to majordomo_at_[hidden], with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/