Julian Cummings wrote:
> Just a couple minor comments about the point below. First, in your example of
> the simplified stencil definition using the pre-defined Laplacian3D function,
> don't you still have to include offsets for the other terms? That is,
> shouldn't your example read
>
> DECLARE_STENCIL4(acoustic3D_stencil, P1, P2, P3, c)
>
> P3(0,0,0) = 2 * P2(0,0,0) + c(0,0,0) * Laplacian3D(P2) - P1(0,0,0);
>
> END_STENCIL
I overloaded some operators so that "P3" is equivalent to P3(0,0,0)
and "P2" is P2(0,0,0). This simplifies the notation somewhat.
> My other comment is that in my view it would be nice if stencil objects were
> truly composable from other stencil objects. Your example differs from that
> ideal because Laplacian3D is an inline function, not a stencil itself. I'd
> like to see something like
>
> DECLARE_STENCIL2(Laplacian3D_stencil, P1, P2)
>
> P2(0,0,0) = -6 * P1(0,0,0) + P1(-1,0,0) + P1(1,0,0) + P1(0,-1,0)
> + P1(0,1,0) + P1(0,0,-1) + P1(0,0,1);
> END_STENCIL
>
> DECLARE_STENCIL5(acoustic3D_stencil, P1, P2, P3, P4, c)
>
> P4(0,0,0) = 2 * P2(0,0,0) + c(0,0,0) * Laplacian3D_stencil(P2,P3)
> - P1(0,0,0);
>
> END_STENCIL
>
> I'm not sure this can really be made to work; just seems like it would be
> nice if it could.
In the current design, stencils contain the assignment; stencil
operators like Laplacian3D have to return a scalar value. If you
wanted to compose actual stencils (including assignment), you could do
something like:
DECLARE_STENCIL5(acoustic3D_stencil, P1, P2, P3, P4, c)
Laplacian3D_stencil::apply(P4,P2);
P3 = 2 * P2 + c(0,0,0) * P4 - P1;
END_STENCIL
This requires a temporary array, though. I think that Laplacian3D
really needs to be a stencil operator, rather than a stencil.
I could provide another macro to declare stencil operators, e.g.
DECLARE_STENCIL_OPERATOR1(Laplacian3D, A)
return -6 * A(0,0,0) + A(-1,0,0) + A(1,0,0) + A(0,-1,0)
+ A(0,1,0) + A(0,0,-1) + A(0,0,1);
END_STENCIL_OPERATOR
which would then allow you to write:
DECLARE_STENCIL4(acoustic3D_stencil, P1, P2, P3, c)
P3 = 2 * P2 + c * Laplacian3D(P2) - P1;
END_STENCIL
Would that be better?
Cheers,
Todd
--------------------- blitz-dev list --------------------------------
* To subscribe/unsubscribe: mail to majordomo@oonumerics.org, with
"subscribe blitz-dev" or "unsubscribe blitz-dev" in the body of the message
* Blitz++ web page: http://oonumerics.org/blitz/
This archive was generated by hypermail 2b29 : Wed Feb 20 2002 - 04:30:05 EST