POOMA: A C++ Toolkit for High-Performance Parallel Scientific Computing | ||
---|---|---|
Prev | Chapter 3. A Tutorial Introduction | Next |
POOMA Arrays support many scientific computations, but other scientific computations require values distributed throughout space, and Arrays have no spatial extent. POOMA Fields, supporting a superset of Array functionality, model values distributed throughout space.
A Field consists of a set of cells distributed through space. Like an Array cell, each Field cell is addressed via indices. Unlike an Array cell, each Field cell can hold multiple values. Like Arrays, Fields can be accessed via data-parallel expressions and stencils and may be distributed across processors. Unlike Array cells, Field cells exist in a multidimensional volume so, e.g., distances between cells and normals to cells can be computed.
In this section, we implement the Doof2d two-dimensional diffusion simulation program using Fields. This simulation does not require any Field-specific features, but we present this program rather than one using Field-specific features to facilitate comparison with the Array versions, especially Example 3-3.
Example 3-6. Data-Parallel Field Implementation of Doof2d
#include <iostream> // has std::cout, ... #include <stdlib.h> // has EXIT_SUCCESS #include "Pooma/Fields.h" // has POOMA's Field declarations// Doof2d: POOMA Fields, data-parallel implementation int main(int argc, char *argv[]) { // Prepare the POOMA library for execution. Pooma::initialize(argc,argv); // Ask the user for the number of averagings. long nuAveragings, nuIterations; std::cout< <"Please enter the number of averagings: "; std::cin > > nuAveragings; nuIterations = (nuAveragings+1)/2; // Each iteration performs two averagings. // Ask the user for the number n of values along // one dimension of the grid. long n; std::cout < < "Please enter the field size: "; std::cin > > n; // Specify the fields' domains [0,n) x [0,n). Interval<1> N(0, n-1); Interval<2> vertDomain(N, N); // Set up interior domains [1,n-1) x [1,n-1) for // computation. Interval<1> I(1,n-2); Interval<1> J(1,n-2); // Specify the fields' mesh, i.e., its spatial // extent, and its centering type.
DomainLayout<2> layout(vertDomain); UniformRectilinearMesh<2> mesh(layout, Vector<2>(0.0), Vector<2>(1.0, 1.0)); Centering<2> cell = canonicalCentering<2>(CellType, Continuous, AllDim); // Create the fields. // The Field template parameters indicate a mesh, a // 'double' value type, and ordinary 'Brick' // storage.
Field<UniformRectilinearMesh<2>, double, Brick> a(cell, layout, mesh); Field<UniformRectilinearMesh<2>, double, Brick> b(cell, layout, mesh); // Set up the initial conditions. // All grid values should be zero except for the // central value. a = b = 0.0; // Ensure all data-parallel computation finishes // before accessing a value. Pooma::blockAndEvaluate(); b(n/2,n/2) = 1000.0; // In the average, weight elements with this value. const double weight = 1.0/9.0; // Perform the simulation. for (int k = 0; k < nuIterations; ++k) { // Read from b. Write to a.
a(I,J) = weight * (b(I+1,J+1) + b(I+1,J ) + b(I+1,J-1) + b(I ,J+1) + b(I ,J ) + b(I ,J-1) + b(I-1,J+1) + b(I-1,J ) + b(I-1,J-1)); // Read from a. Write to b. b(I,J) = weight * (a(I+1,J+1) + a(I+1,J ) + a(I+1,J-1) + a(I ,J+1) + a(I ,J ) + a(I ,J-1) + a(I-1,J+1) + a(I-1,J ) + a(I-1,J-1)); } // Print out the final central value. Pooma::blockAndEvaluate(); // Ensure all computation has finished. std::cout < < (nuAveragings % 2 ? a(n/2,n/2) : b(n/2,n/2)) < < std::endl; // The fields are automatically deallocated. // Tell the POOMA library execution has finished. Pooma::finalize(); return EXIT_SUCCESS; }
As mentioned above, the fundamental difference between Arrays and Fields is the latter has cells and meshes. The Field declarations reflect this. To declare a Field, the Pooma/Fields.h header file must be included. A Field's domain consists of a set of cells, sometimes called positions when referring to Arrays. As for Arrays, a Field's domain and its layout must be specified. Since the above program is designed for uniprocessor computation, specifying the domain specifies the layout. A Field's mesh specifies its spatial extent. For example, one can ask the mesh for the distance between two cells or for the normals to a particular cell. Cells in a UniformRectilinearMesh all have the same size and are parallelepipeds. To create the mesh, one specifies the layout, the location of the spatial point corresponding to the lower, left domain location, and the size of a particular cell. Since this program does not use mesh computations, our choices do not matter. We specify the domain's lower, left corner as spatial location (0.0, 0.0) and each cell's width and height as 1. Thus, the middle of the cell at domain position (3,4) is (3.5, 4.5).
A Field cell can contain one or more values although each cell must have the same arrangement of values. For this simulation, we desire one value per cell so we place that position at the cell's center, i.e., a cell centering. The canonicalCentering function returns such a centering. .
A Field declaration is analogous to an Array declaration but must also specify a centering and a mesh. In Example 3-3, the Array declaration specifies the array's dimensionality, the value type, the Engine type, and a layout. Field declarations specify the same values. Its first template parameter specifies the mesh's type, which includes an indication of its dimensionality. The second and third template parameters specify the value type and the Engine type. Since a Field has a centering and a mesh in addition to a layout, those arguments are also necessary.
Field operations are a superset of Array operations so the Doof2d computations are the same as in Example 3-3. Field accesses require parentheses, not square brackets, and accesses to individual values should be preceded by calls to Pooma::blockAndEvaluate.
To summarize, Fields support multiple values per cell and have spatial extent. Thus, their declarations must specify a centering and a mesh. Otherwise, a Field program is similar to one using Arrays.