|
|
Blitz Support : |
From: Andreas R. (andreasreifschneider_at_[hidden])
Date: 2005-05-31 07:17:33
Faheem Mitha wrote:
>
>
> On Mon, 30 May 2005, Faheem Mitha wrote:
>
>> Also, I was very disappointed to discover that stencils do not accept
>
> scalar
>
>> values. This is a *major* usability problem. At least for me, this
>> means I won't be able to use stencils, since working around this in my
>> context would be so difficult as to be next to impossible.
>>
>> I could not find much discussion about this, but in examples/cfd.cpp,
>> I read:
>>
>> * - Stencil objects can't take scalar arguments, so all the scalars
>> * have to be globals. Big ugh. There is a fix for this,
>> * but it hasn't been implemented yet.
>>
>> Apparently it still has not been implemented, since I get errors like
>> " error: `const double' is not a class, struct, or union type" when I
>> try to pass values as an argument.
>>
>> Are there any plans to fix this, and what was the fix referred to?
>
>
> I came across what appears to be a workaround,
> http://www.oonumerics.org/MailArchives/blitz-support/2004/06/1128.php
>
> Is this the best one can do?
If you have many stencils you can factor out the common code into a base
class, e.g. like this:
/** @param boundary_size
It is a template parameter, since this is faster than
passing it through constructor.
Apart from that this wouldn't bring any flexibility,
since the Blitz++ stencil operators require compile time specification
of it (one has to chosse whether one uses Laplacian2D or Laplacian2D4).
If @a boundary_size == 1 one has to use the 2nd order accurate
approximation in the stencil (e.g. Laplacian2D);
If @a boundary_size == 2 one has to use the 4th order accurate
approximation in the stencil (e.g. Laplacian2D4).
@todo generalize the Stencils to accept boundary_size
*/
template <int boundary_size = 1>
struct Stencil {
// Stencil(int boundary_size = 1) : boundary_size( boundary_size ) {}
enum { BOUNDARY_SIZE = boundary_size };
template <int N>
inline void getExtent(TinyVector<int,N>& minb,
TinyVector<int,N>& maxb) const {
minb = shape(-boundary_size,-boundary_size);
maxb = shape(+boundary_size,+boundary_size);
}
enum { hasExtent = 1 };
// private:
// int boundary_size;
};
The actual stencil is then somewhat shorter (just an example from my code):
struct SimpleStripes2D_a_Stencil : public Stencil<1> {
SimpleStripes2D_a_Stencil(double d_a, double b_a, double r_a, double s_a)
: d_a(d_a), b_a(b_a), r_a(r_a), s_a(s_a) {}
template <typename T1, typename T2, typename T3, typename T4,
typename T5, typename T6, typename T7, typename T8,
typename T9, typename T10, typename T11>
inline void apply(T1& a_out, T2& a, T3& b, T4& s_beta, T5&,
T6&, T7&, T8&, T9&, T10&, T11&) const
{
a_out = a + s_beta * (a*a + b_a)/(b * (1 + s_a*a*a)) - r_a*a +d_a *
Laplacian2D(a);
if (a_out < 0.0) a_out = 0.0;
}
private:
const double d_a;
const double b_a, r_a, s_a;
};
Regards,
Andreas