Proto
Data Structures | Public Member Functions
Proto::Stencil< T > Class Template Reference

A Linear Stencil Operation. More...

#include <Proto_Stencil.H>

Data Structures

struct  coeff_holder
 

Public Member Functions

template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
void apply (const BoxData< T, C, MEMTYPE, D, E > &a_src, BoxData< T, C, MEMTYPE, D, E > &a_dest, const Box &a_bx, bool a_initToZero=false, const T a_scale=1) const
 Apply Stencil Helper function. More...
 
Constructors
 Stencil ()
 Default Constructor.
 
 Stencil (Shift a_shift, T a_coef, Point a_destRefratio=Point::Ones(), Point a_destShift=Point::Zeros(), Point a_srcRefratio=Point::Ones())
 General Constructor. More...
 
Operators
Stencil< T > operator* (const Stencil< T > &a_stencil) const
 Stencil Composition. More...
 
Stencil< T > operator* (const T a_coef) const
 Scalar Multiplication. More...
 
void operator*= (const Stencil< T > &a_stencil)
 In Place Stencil Composition. More...
 
void operator*= (const T a_coef)
 In Place Scalar Multiplication.
 
Stencil< T > operator+ (const Stencil< T > &a_stencil) const
 Stencil Addition. More...
 
Stencil< T > operator- (const Stencil< T > &a_stencil) const
 Stencil Subtraction (Convenience) More...
 
void operator+= (const Stencil< T > &a_stencil)
 In Place Stencil Addition. More...
 
void operator-= (const Stencil< T > &a_stencil)
 In Place Stencil Subtraction. More...
 
Stencil Binding and Application
template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
LazyStencil< T, C, MEMTYPE, D, E > operator() (const BoxData< T, C, MEMTYPE, D, E > &a_src, T a_scale=1) const
 Operate on BoxData. More...
 
template<unsigned int C, MemType MEMTYPE, unsigned char D, unsigned char E>
LazyStencil< T, C, MEMTYPE, D, E > operator() (const BoxData< T, C, MEMTYPE, D, E > &a_src, Box a_box, T a_scale=1) const
 Operate on BoxData (Overload with Box Input) More...
 
Utility
void invert (int a_dir)
 Invert Stencil. More...
 
void transpose (unsigned char a, unsigned char b)
 Transpose Stencil. More...
 
Box indexRange (Box a_domain) const
 Get Max Index Range. More...
 
Box indexDomain (Box a_range) const
 Get Max Index Domain. More...
 
Box range (Box a_domain) const
 Get Max Range Box. More...
 
Box domain (Box a_range) const
 Get Min Domain Box. More...
 
diagonalValue () const
 Get Diagonal Value. More...
 
void print () const
 Print. More...
 

Static Public Member Functions

Stencil Library
static Stencil< T > Derivative (int a_n, int a_dir=0, int a_order=2)
 Stencil Library: Derivative. More...
 
static Stencil< T > Laplacian ()
 Stencil Library: Laplacian. More...
 
static Stencil< T > LaplacianFace (int a_dir, int a_order=2)
 Stencil Library: Perpendicular Laplacian. More...
 
static Stencil< T > CellToEdge (int a_dir, int a_order=4)
 Stencil Library: Cell to Edge Deconvolution. More...
 
static Stencil< T > DiffCellToEdge (int a_dir, int a_order=4)
 Stencil Library: Cell to Edge Differentiation. More...
 
static Stencil< T > EdgeToCell (int a_dir, int a_order=4)
 Stencil Library: Edge to Cell Convolution. More...
 
static Stencil< T > CellToEdgeL (int a_dir, int a_order=5)
 Stencil Library: Downwind Cell to Edge Deconvolution. More...
 
static Stencil< T > CellToEdgeH (int a_dir, int a_order=5)
 Stencil Library: Upwind Cell to Edge Deconvolution. More...
 
static Stencil< T > AvgDown (int a_refRatio=2)
 Stencil Library: Simple Average. More...
 
static Stencil< T > AvgDownEdge (int a_dir, int a_refRatio=2)
 Stencil Library: Simple Average over an Edge. More...
 
static Stencil< T > FluxDivergence (int a_dir)
 Stencil Library: Flux Divergence. More...
 

Accessors and Queries

bool operator== (Stencil< T > &a_stencil) const
 
bool operator!= (Stencil< T > &a_stencil) const
 Inquality Operator. More...
 
const std::vector< T > & coefs () const
 Get Vector of Coefficients. More...
 
const T * devCoefs () const
 Get Vector of Coefficients. More...
 
const PointdevOffsets () const
 Get Vector of Offsets. More...
 
const std::vector< Point > & offsets () const
 Get Vector of Offsets. More...
 
unsigned long long int size () const
 Size. More...
 
Box span () const
 Span. More...
 
Point spanPoint () const
 Span. More...
 
Point ghost () const
 Ghost. More...
 
PointsrcRatio ()
 Get Source Refinement Ratio. More...
 
const PointsrcRatio () const
 Get Source Refinement Ratio (Const)
 
PointdestRatio ()
 Get Destination Refinement Ratio. More...
 
const PointdestRatio () const
 Get Destination Refinement Ratio (Const)
 
PointdestShift ()
 Get Destination Shift. More...
 
const PointdestShift () const
 Get Destination Shift (Const)
 

Detailed Description

template<typename T>
class Proto::Stencil< T >

A Linear Stencil Operation.

Encapsulates a linear stencil operation where coefficients are of type T. Stencil objects are built and used in a way that conforms to their nature as operators. For illustrative usage examples, refer to the following code snippets:

Examples: Build a Stencil from Shifts and coefficients:

// This example builds and applies the standard 2*DIM + 1 Point Laplacian Stencil:
// In 1D, this can be done in a single line:
// stencil = 1 x low side - 2 x center + 1 x high side
Stencil<double> L_1D = 1.0*Shift::Basis(0,-1) - 2.0*Shift::Zeros() + 1.0*Shift::Basis(0,1);
// In arbitrary dimension:
Stencil<double> L = (-2.0*DIM)*Shift::Zeros();
for (int dir = 0; dir < DIM; dir++)
{
L += 1.0*Shift::Basis(dir,1) + 1.0*Shift::Basis(dir,-1);
}

Apply a Stencil with no Source / Dest Refinement to a BoxData:

// Build a Stencil from the built in Proto Stencil library:
auto L = Stencil<double>::Laplacian(); //2nd order 2*DIM + 1 Point Laplacian
// Initialize some source data
// For this example, we define our input based on our desired ouput: an 8^DIM Point cube.
int rangeSize = 8;
Box rangeBox = Box::Cube(rangeSize); // [(0,...,0), (7,...,7)]
// The domainBox can be computed from the desired range, tanking ghost (halo) cells into account
Box domainBox = L.domain(rangeBox); // L needs 1 layer of ghost cells, so the domainBox is: [ (-1,...,-1), (8,....,8)]
Box computeBox = rangeBox.grow(-1);
// Initialize the data using forall
double dx = 0.5*M_PI/rangeSize; //assuming a range of PI/2
BoxData<double> Src = forall_p<double>([=] PROTO_LAMBDA (Point& pt, Var<double>& src)
{
src(0) = sin(pt[0]*dx); //set Src = sin(x)
}, domainBox); //don't forget the domain Box!
// Apply L to Src
// Method 1: Build a new BoxData from the Stencil computation
BoxData<double> Dest_0 = L(Src);
// When no Box input is specified, the output will be defined on the largest possible Box (e.g. rangeBox)
BoxData<double> Dest_1 = L(Src, computeBox);
// When a VALID Box input is specified, the output will be restricted to that Box.
// *Usually* it makes more sense to let Stencil do its magic and range Box for you.
// Method 2: Apply Stencil computation output to an existing BoxData
BoxData<double> Dest_3(rangeBox);
BoxData<double> Dest_4(rangeBox);
Dest_4.setVal(0);
// REPLACE the data in the destination with the output of L(Src)
Dest_3 |= L(Src);
// ADD the output of L(Src) to the existing data in the destination
Dest_4 += L(Src);
// Both the ADD and REPLACE operations can specify a compute Box as well if desired
// Again, specifying the Box input is not recommended unless it is necessary
// WARNING: Note that 'auto' does NOT work for Stencil evaluation!!
// Compile errors involving 'LazyStencil' could be caused by inappropriate use of 'auto'.

The above examples illustrate the usage of Stencil to do computations in which the source and destination arrays are at the same refinement. Stencil is also capable of "prolong" and "restrict" type operations in which the destination and source respectively are refined. This functionality is useful for operations such as averaging and interpolation, or for algorithms like Red-Black Gauss Seidel Iteration in which it is necessary to evaluate a Stencil on "every other cell".

To facilitate these more exotic operations, the Stencil API allows the user to designate a source and/or destination refinement ratio. If these values are different from (1,...,1), then input Box object will be interpreted as an "index range" instead of a physical Box domain. The following code snippets illustrate some examples of this functionality.

Examples: Non-Trivial Source Refinement Ratio

// This example Stencil computes a linear average from fine data source data onto a coarse grid
Stencil<double> Avg;
int refRatio = 2;
Box offsetBox = Box::Cube(refRatio);
for (auto iter = offsetBox.begin(); iter != offsetBox.end(); ++iter)
{
Avg += 1.0*Shift(*iter);
}
Avg *= (1.0/offsetBox.size());
// WARNING: Stencils with different src/dest refinement ratios cannot be added, subtracted, etc.
// When building a Stencil with non-trivial refinement, set the refinement ratio last.
Avg.srcRatio() = Point::Ones(refRatio);
int rangeSize = 8;
Box indexBox = Box::Cube(rangeSize); //[(0,...,0), (7,....,7)]
Box domainBox = indexBox.refine(refRatio); //[(0,...,0), (15,...,15)]
auto Src = forall_p<double>([] PROTO_LAMBDA (Point& pt, Var<double>& src)
{
src(0) = 0.0;
for (int ii = 0; ii < DIM; ii++)
{
src(0) += pt[ii];
}
},domainBox);
// The index box represents the points at which to compute the Stencil.
// When the source is refined by r, the data used to compute the solution at point p and p + 1
// In the index box is shifted by r.
// In this example, r = 2. In 1D, the average at Point (0) is computed by:
// avg(0) = (1.0*src(0) + 1.0*src(1))/2.
// At the next Point, the forumula is:
// avt(1) = (1.0*src(2) + 1.0*src(3))/2.
// Manifesting a shift in the source data by r = 2.
BoxData<double> Dest_0 = Avg(Src,indexBox);
// OR
BoxData<double> Dest_1 = Avg(Src); //Stencil automatically determines the largest possible Box for Dest_1, given the data available.
// The "|=" and "+=" operators can be used here as well with identical semantics.
// For an example illustrating the usage of the destination refinement ratio, see the documentation for InterpStencil

Non-Trivial Destination Refinement Ratio

int refRatio = 2;
// This is a very simple stencil; it scales a Point by 7.
Stencil<double> S = 7.0*Shift::Zeros();
S.destRatio() = Point::Ones(2);
S.destShift() = Point::Ones();
int domainSize = 8;
Box indexBox = Box::Cube(domainSize); //[(0,...,0), (7,....,7)]
Box rangeBox = indexBox.refine(refRatio); //[(0,...,0), (15,...,15)]
BoxData<double> Src(indexBox);
Src.setVal(1);
BoxData<double> Dest(rangeBox);
Dest.setVal(0);
Dest |= S(Src);
// In the case of a non-trivial source refinement ratio, the Stencil is shifted by srcRatio instead of by 1 for each
// consecutive Stencil evaluation. By contrast, when the destination refinement ratio is non-trivial, the Stencil only
// populates one in every destRefRatio Points, with the remainder of points being untouched by the operation.
// In this particular example, at each Point in the indexBox, the following graphic represents the Stencil update (DIM = 2):
//
// SOURCE DESTINATION DESTINATION
// +-------------+ +------+------+ +------+------+
// | | | | | | | |
// | | | 0 | 0 | | 0 | 7 |
// | 1 | S +------+------+ ---> +------+------+
// | | | | | | | |
// | | | 0 | 0 | | 0 | 0 |
// +-------------+ +------+------+ +------+------+
//

In the case of non-trivial destination refinement, an array of refRatio^DIM Stencils is often needed to fully populate the output. Proto provides a convenience structure called InterpStencil which is designed to mitigate this form of pedantry. See the associated documentation for additional information and example usage.


The documentation for this class was generated from the following files: