Index: StatigraphicFlow.cpp =================================================================== RCS file: /home/pooma/Repository/r2/examples/NewField/StatigraphicFlow/Attic/StatigraphicFlow.cpp,v retrieving revision 1.1.2.3 diff -c -p -r1.1.2.3 StatigraphicFlow.cpp *** StatigraphicFlow.cpp 2001/08/23 23:01:20 1.1.2.3 --- StatigraphicFlow.cpp 2001/08/29 16:03:28 *************** *** 1,5 **** // Oldham, Jeffrey D. ! // 2001Jul25 // Pooma // Chevron Kernel Written Using POOMA's Proposed Field Abstraction --- 1,5 ---- // Oldham, Jeffrey D. ! // 2001Aug23 // Pooma // Chevron Kernel Written Using POOMA's Proposed Field Abstraction *************** *** 107,129 **** /** THE PROGRAM **/ ! // Return the index of the specified field offset in the given list. ! // Return a negative number if not found. template ! inline int ! findIndex(const FieldOffsetList &vec, ! const FieldOffset &fo) { ! int indx; ! for (indx = vec.size()-1; ! indx >= 0 && vec[indx] != fo; ! --indx) ! ; ! return indx; ! } int main(int argc, char *argv[]) { // Set up the Pooma library. --- 107,287 ---- /** THE PROGRAM **/ ! // Use a linear algebra computation that stores the pressure gradient ! // values in the four quadrants around the current vertex (loc). template ! struct ComputeGradientsInfo { ! void scalarCodeInfo(ScalarCodeInfo &info) const ! { ! // ComputeGradients::operator() has 4 arguments: vertices to ! // iterate over and three fields. ! ! info.arguments(5); ! ! // Only the first such argument is changed. ! ! info.write(0, false); ! info.write(1, true); ! info.write(2, false); ! info.write(3, false); ! info.write(4, false); ! ! // Does the operation index neighboring cells, i.e., do we need to ! // update the internal guard layers? ! ! info.useGuards(0, false); ! info.useGuards(1, true); ! info.useGuards(2, true); ! info.useGuards(3, true); ! info.useGuards(4, true); ! ! info.dimensions(Dim); ! ! for (int i = 0; i < Dim; ++i) { ! info.lowerExtent(i) = 1; // We access neighbors. ! info.upperExtent(i) = 0; ! } ! } ! }; ! ! ! template ! struct ComputeGradients ! : public ComputeGradientsInfo ! { ! ! ComputeGradients(const Centering &disspoke, ! const FieldOffsetList &gradients, ! const int nuFluxPoints, ! const std::vector > &disFluxPoints, ! const Centering &subcell) ! : disspoke_m(disspoke), gradients_m(gradients), ! nuFluxPoints_m(nuFluxPoints), disFluxPoints_m(disFluxPoints), ! subcell_m(subcell) ! { ! PAssert(gradients_m.size() * Dim == nuFluxPoints_m * 2); ! } ! ! // FIXME: Perhaps we want to modify ScalarCode to take a first ! // FIXME: argument of a centering. In the meantime, we just use a ! // FIXME: fake field. ! ! template ! inline ! void operator()(const F1 &vertexField, ! F2 &pressureGradient, ! const F3 &faceDistance, ! const F4 &directedPermeability, ! const F5 &pressure, ! const Loc &loc) const ! { ! const int nuRows = (1 << Dim) * Dim; ! TNT::Matrix A(nuRows, nuRows, 0.0); ! TNT::Vector rhs(nuRows, 0.0); ! TNT::Vector ipiv; // ignored ! ! // Assign values to the matrix A and vector rhs. ! ! for (int faceIndex = nuFluxPoints_m-1; faceIndex >= 0; --faceIndex) { ! // Work on the "positive" side of the face. ! FieldOffset fo = disFluxPoints_m[faceIndex][0]; ! int columnNu = ! findIndex(gradients_m, ! nearestNeighbors(pressureGradient.centering(), ! fo, disspoke_m, true)[0]); ! // The column number is the pressure gradient corresponding to the ! // "positive" side of the face. ! ! for (int i = 0; i < Dim; ++i) ! A[faceIndex][columnNu] = ! directedPermeability(nearestNeighbors(directedPermeability.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset())(i); ! A[faceIndex+nuFluxPoints_m][columnNu] = ! faceDistance(nearestNeighbors(faceDistance.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset()); ! rhs[faceIndex+nuFluxPoints_m] -= ! pressure(nearestNeighbors(pressure.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset()); ! ! // Work on the "negative" side of the face. ! fo = disFluxPoints_m[faceIndex][1]; ! columnNu = ! findIndex(gradients_m, ! nearestNeighbors(pressureGradient.centering(), ! fo, disspoke_m, true)[0]); ! // The column number is the pressure gradient corresponding to the ! // "positive" side of the face. ! ! for (int i = 0; i < Dim; ++i) ! A[faceIndex][columnNu] = ! -directedPermeability(nearestNeighbors(directedPermeability.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset())(i); ! A[faceIndex+nuFluxPoints_m][columnNu] = ! -faceDistance(nearestNeighbors(faceDistance.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset()); ! rhs[faceIndex+nuFluxPoints_m] -= ! -pressure(nearestNeighbors(pressure.centering(), ! fo, disspoke_m, true)[0], ! loc+fo.cellOffset()); ! } ! ! // Solve for the pressure gradients. ! ! TNT::LU_solve(A, ipiv, rhs); ! ! // Now, rhs has the pressure gradients. ! ! for (int faceIndex = nuFluxPoints_m-1; faceIndex >= 0; --faceIndex) ! for (int i = 0; i < Dim; ++i) ! pressureGradient(gradients_m[faceIndex], loc)(i) = rhs[faceIndex+i]; ! ! return; ! } ! ! // Return the index of the specified field offset in the given list. ! // Return a negative number if not found. ! ! static ! inline int ! findIndex(const FieldOffsetList &vec, ! const FieldOffset &fo) ! { ! int indx; ! for (indx = vec.size()-1; ! indx >= 0 && vec[indx] != fo; ! --indx) ! ; ! return indx; ! } ! ! private: + // Discontinuous spokes. + const Centering &disspoke_m; + // The pressure gradients. + const FieldOffsetList &gradients_m; + + // The number of flux points for a cell. + const int nuFluxPoints_m; + + // For every i, disFluxPoints_m[i] is two discontinuous positions on + // "either side" of the face represented by the flux points + // fluxPoints[i]. + const std::vector > &disFluxPoints_m; + + // One cell value in each quadrant. + const Centering &subcell_m; + }; + + int main(int argc, char *argv[]) { // Set up the Pooma library. *************** int main(int argc, char *argv[]) *** 172,178 **** position(1) = 0.75; subcell.addValue(Orientation(1), position); position(0) = 0.75; subcell.addValue(Orientation(1), position); position(1) = 0.25; subcell.addValue(Orientation(1), position); ! Fields_t pressureGradient (subcell, meshLayout, origin, spacings); // Spoke-centered Field. --- 330,336 ---- position(1) = 0.75; subcell.addValue(Orientation(1), position); position(0) = 0.75; subcell.addValue(Orientation(1), position); position(1) = 0.25; subcell.addValue(Orientation(1), position); ! Fieldv_t pressureGradient (subcell, meshLayout, origin, spacings); // Spoke-centered Field. *************** int main(int argc, char *argv[]) *** 203,211 **** Centering disFace = canonicalCentering(FaceType, Discontinuous); Fieldv_t directedPermeability (disFace, meshLayout, origin, spacings); ! // \gamma_{i,j} = K_i^t \dot \hat{n}_j Fields_t faceDistance (disFace, meshLayout, origin, spacings); ! // distance from cell center to face center /* INITIALIZATION */ --- 361,369 ---- Centering disFace = canonicalCentering(FaceType, Discontinuous); Fieldv_t directedPermeability (disFace, meshLayout, origin, spacings); ! // \gamma_{i,j} = K_i^t \dot \hat{n}_j Fields_t faceDistance (disFace, meshLayout, origin, spacings); ! // distance from cell center to face center /* INITIALIZATION */ *************** int main(int argc, char *argv[]) *** 238,325 **** faceDistance.centering()); #endif // PSEUDOCODE - const int nuRows = (1 << Dim) * Dim; - TNT::Matrix A(nuRows, nuRows, 0.0); - TNT::Vector rhs(nuRows, 0.0); - TNT::Vector ipiv; // ignored - - // FIXME: Move this code to a stencil so it can be applied across - // the entire grid. - - // Assign values to the matrix A and vector rhs. - const Centering vert = canonicalCentering(VertexType, Continuous); PInsist(vert.size() == 1, "The vertex centering has too many values."); const FieldOffsetList gradients = nearestNeighbors(pressureGradient.centering(), FieldOffset(Loc(0)) /* cell origin */, vert); ! // gradients's order of pressure gradients will be used for the ! // matrix and rhs. const FieldOffsetList fluxPoints = nearestNeighbors(spokeFlux.centering(), FieldOffset(Loc(0)) /* cell origin */, vert); ! // fluxPoints has locations for all faces incident to the vertex. ! const int size = fluxPoints.size(); ! PAssert(gradients.size() * Dim == size * 2); const std::vector > disFluxPoints = nearestNeighbors(disspoke, fluxPoints, spokeFlux.centering()); ! // For every i, disFluxPoints[i] is two discontinuous positions on ! // "either side" of the face represented by fluxPoints[i]. ! ! for (int faceIndex = size-1; faceIndex >= 0; --faceIndex) { ! // Work on the "positive" side of the face. ! FieldOffset fo = disFluxPoints[faceIndex][0]; ! int columnNu = ! findIndex(gradients, ! nearestNeighbors(pressureGradient.centering(), ! fo, disspoke, true)[0]); ! // The column number is the pressure gradient corresponding to the ! // "positive" side of the face. ! ! // FIXME: The lhs (double) and rhs (vector field) do not match. ! ! A[faceIndex][columnNu] = ! directedPermeability(nearestNeighbors(directedPermeability.centering(), ! fo, disspoke, true)[0]); ! A[faceIndex+size][columnNu] = ! faceDistance(nearestNeighbors(faceDistance.centering(), ! fo, disspoke, true)[0]); ! rhs[faceIndex+size] -= ! pressure(nearestNeighbors(pressure.centering(), ! fo, disspoke, true)[0]); ! ! fo = disFluxPoints[faceIndex][1]; ! columnNu = ! findIndex(gradients, ! nearestNeighbors(pressureGradient.centering(), ! fo, disspoke, true)[0]); ! // The column number is the pressure gradient corresponding to the ! // "positive" side of the face. ! ! A[faceIndex][columnNu] = ! -directedPermeability(nearestNeighbors(directedPermeability.centering(), ! fo, disspoke, true)[0]); ! A[faceIndex+size][columnNu] = ! -faceDistance(nearestNeighbors(faceDistance.centering(), ! fo, disspoke, true)[0]); ! rhs[faceIndex+size] -= ! -pressure(nearestNeighbors(pressure.centering(), ! fo, disspoke, true)[0]); ! } ! ! // Solve for the pressure gradients. ! ! TNT::LU_solve(A, ipiv, rhs); ! ! // Now, rhs has the pressure gradients. ! for (int faceIndex = size-1; faceIndex >= 0; --faceIndex) ! // FIXME: Is this type of assignment supported by the current code base? ! pressureGradient(gradients[faceIndex], subcell) = rhs[faceIndex]; // Compute the spoke fluxes. --- 396,433 ---- faceDistance.centering()); #endif // PSEUDOCODE const Centering vert = canonicalCentering(VertexType, Continuous); PInsist(vert.size() == 1, "The vertex centering has too many values."); const FieldOffsetList gradients = nearestNeighbors(pressureGradient.centering(), FieldOffset(Loc(0)) /* cell origin */, vert); ! // gradients's order of pressure gradients will be used for the ! // matrix and rhs. const FieldOffsetList fluxPoints = nearestNeighbors(spokeFlux.centering(), FieldOffset(Loc(0)) /* cell origin */, vert); ! // fluxPoints has locations for all faces incident to the vertex. ! ! const int nuFluxPoints = fluxPoints.size(); const std::vector > disFluxPoints = nearestNeighbors(disspoke, fluxPoints, spokeFlux.centering()); ! // For every i, disFluxPoints[i] is two discontinuous positions on ! // "either side" of the face represented by fluxPoints[i]. ! typedef ComputeGradients CG_t; ! CG_t cG(disspoke, gradients, nuFluxPoints, disFluxPoints, subcell); ! ScalarCode computeGradients(cG); ! ! // FIXME: Use an otherwise unused field for the ScalarCode ! // FIXME: iteration over vertices. ! Fields_t vertexField(canonicalCentering(VertexType, Continuous), ! meshLayout, origin, spacings); + computeGradients(vertexField, pressureGradient, + faceDistance, directedPermeability, pressure); // Compute the spoke fluxes. *************** int main(int argc, char *argv[]) *** 346,359 **** totalFlux = sum(spokeFlux.mesh().normals().signedMagnitude() * - // FIXME: This is not yet implemented. We want a - // FIXME: data-parallel sum. This is we want a function - // FIXME: Field_t sum(/* input */ Field_t, - // FIXME: std::vector, /*output */ - // FIXME: Centering). The vector's length == the output - // centering's length. The function works by using the input - // field with each FieldOffsetList to form one value in the - // output field. sum(spokeFlux, nearestNeighbors(spokeFlux.centering(), disFace), disFace), --- 454,459 ----