[pooma-dev] Chevron Code Using New Field Abstractions

[Jeffrey Oldham]

On Wed, Jul 25, 2001 at 02:12:49PM -0700, Jeffrey Oldham wrote:
> Attached is a very preliminary version of the Chevron code written
> using C++ pseudocode closely related to the proposed NewField
> revisions.  It does not compile since the underlying NewField and mesh
> routines have not yet been implemented.
> 
> The next steps are:
> 
> 1. To ensure that the algorithm is correct.
> 2. To add more comments describing my assumptions about functions and
>    classes.

I have added comments near the beginning of the file.

> 3. To discuss whether the syntax is acceptable.
> 4. To make the code available in some portion of the Pooma CVS tree.

Thanks,
Jeffrey D. Oldham
oldham at codesourcery.com
-------------- next part --------------
// Oldham, Jeffrey D.
// 2001Jul25
// Pooma

// Chevron Kernel Written Using POOMA's Proposed Field Abstraction

#include <iostream>
#include <cstdlib>
#include "Pooma/NewFields.h"

// This program implements "Implementation of a Flux-Continuous Fnite
// Difference Method for Stratigraphic, Hexahedron Grids," by
// S. H. Lee, H. Tchelepi, and L. J. DeChant, \emph{1999 SPE Reservoir
// Simulation Symposium}, SPE (Society of Petroleum Engineers) 51901.

// Preprocessor symbols:
// PSEUDOCODE: Do not define this symbol.  Surrounds desired code to
//	       deal with different granularity fields.
// DEBUG: If defined, print some information about internal program
//	  values.


/** QUESTIONS **/

// o. If several different fields are created using the same mesh
//    object, is the mesh object shared?
// o. Can meshes be queried without going through an associated field?
// o. According to my understanding, the Chevron algorithm should be
//    imbedded inside a loop of some type that repeatedly updates the
//    coordinates.
// o. I omitted a separate coordinates field, presumably updated each
//    iteration, in favor of using the mesh.  Since I do not know how
//    the coordinates are updated, I omitted updating the mesh.
// o. Is it important to flesh out the linear algebra solution?  We
//    might learn something about field syntax, but it will also take
//    time for me to determine the correct operands.
// o. The eight spoke-centered flux values are discontinuous, right?
// o. Creating non-canonical edge and face centerings requires
//    dimension-dependent code.  Is this acceptable?


/** UNFINISHED WORK **/

// o ConstField = a Field with values that do not change
// o nearestNeighbors(inputCentering, outputCentering)
// o replicate(field, std::vector<FieldOffsetList>)
// o meshLayout.unitCoordinateNormals()
// o field.mesh()
// o field.mesh().normals()
// o field.mesh().normals().signedMagnitude()
// o sum(field, FieldOffsetList)


/** EXPLANATIONS **/

// o Centering<Dim> canonicalCentering<Dim>(CellType, Continuous):
//    returns a centering object for a cell-centered field with one
//    value at the cell's center (in logical coordinate space)
// o subcell: This centering contains four cell-centered values at
//    positions (0.25, 0.25), (0.25, 0.75), (0.75, 0.75), (0.75, 0.25).  
//    Since this centering is not a canonical centering, it must be
//    constructed.  To do so, we start with a cell-centered centering
//    without any values and repeatedly add values.  The orientation,
//    ignored for cell-centered values, indicates which coordinate values
//    are fixed and which are not.  Using a (1,...,1) indicates that
//    all coordinate values may be changed.
// o spoke: This face-centering has two values on each face.  It, too,
//    has to be constructed since it is not a normal centering.
// o The Chevron algorithm first solves a linear program.  I have
//    omitted since computation since it does not illustrate field
//    computations.
// o replicate(field, std::vector<FieldOffsetList>): This function,
//    syntactic sugar for a nearest neighbors computation, copies the
//    field values to the positions indicated by the
//    std::vector<FieldOffsetList>.  Each field value is copied to one
//    or more values.  replicate() could be replaced by sum(), but the
//    latter function has an unnecessary loop since each output value
//    equals one input value.
// o nearestNeighbors(inputCentering, outputCentering): This function
//    returns a std::vector of FieldOffsetList's, one for each output
//    value specified by the given output centering.  For each output
//    value, the closest input values, wrt Manhattan distance, are
//    returned.  Eventually, these may be pre-computed or cached to
//    reduce running time.
// o meshLayout.unitCoordinateNormals(): This returns a discontinuous
//    face-centered field with unit-length normals all pointing in
//    positive directions.
// o field.mesh(): Returns the mesh object associated with the field.
// o spokeFlux.mesh().normals(): Returns a face-centered field of
//    normal vectors perpendicular to each face.  The magnitude of each
//    normal equals the face's area/volume.
// o spokeFlux.mesh().normals().signedMagnitude(): Returns a
//    face-centered field of scalars, each having absolute value
//    equalling the face's area/volume and sign equalling whether the
//    face's normal is in a positive direction, e.g., the positive
//    x-direction vs. the negative x-direction.
// o sum(field, FieldOffsetList): this parallel-data statement adds
//    the values indicated in the FieldOffsetList to form each output value


/** THE PROGRAM **/

int main(int argc, char *argv[])
{
  // Set up the Pooma library.
  Pooma::initialize(argc,argv);
#ifdef DEBUG
  std::cout << "Start program." << std::endl;
#endif // DEBUG


  /* DECLARATIONS */

  // Create a simple layout.
  const unsigned Dim = 2;		// Work in a 2D world.
  const unsigned nXs = 5;		// number of horizontal vertices
  const unsigned nYs = 4;		// number of vertical vertices
  Interval<Dim> meshDomain;
  meshDomain[0] = Interval<1>(nXs);
  meshDomain[1] = Interval<1>(nYs);
  DomainLayout<Dim> meshLayout(meshDomain, GuardLayers<Dim>(1));

  // Preparation for Field creation.

  Vector<Dim> origin(0.0);
  Vector<Dim> spacings(1.0,1.0);
  typedef UniformRectilinear<Dim, double, Cartesian<Dim> > Geometry_t;
  typedef Field<Geometry_t, double, Brick> Fields_t;
  typedef Field<Geometry_t, double, Brick> ConstFields_t; // TODO: Change to ConstField when ConstField is available.
  typedef Tensor<Dim,double,Full> Tensor_t;
  typedef Field<Geometry_t, Tensor_t, Brick> FieldT_t;
  typedef Field<Geometry_t, Tensor_t, Brick> ConstFieldT_t; // TODO: Change to ConstField when ConstField is available.
  typedef Field<Geometry_t, Vector<Dim>, Brick> Fieldv_t;
  typedef Field<Geometry_t, Vector<Dim>, Brick> ConstFieldv_t; // TODO: Change to ConstField when ConstField is available.

  // Cell-centered Fields.

  Centering<Dim> cell = canonicalCentering<Dim>(CellType, Continuous);
  ConstFieldT_t permeability	(cell, meshLayout, origin, spacings);
  ConstFields_t pressure	(cell, meshLayout, origin, spacings);
  Fields_t totalFlux		(cell, meshLayout, origin, spacings);

  // Subcell-centered Field.

  typedef Centering<Dim>::Orientation Orientation;
  typedef Centering<Dim>::Position Position;
  Position position;
  Centering<Dim> subcell(CellType, Continuous);
  position(0) = 0.25;
  position(1) = 0.25; subcell.addValue(Orientation(1), position);
  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.

  Centering<Dim> spoke(FaceType, Discontinuous);
  // QUESTION: These are supposed to be Discontinuous, right?
  Orientation orientation;
  // NOTE: This code is not dimension-independent.
  for (int zeroFace = 0; zeroFace < 2; ++zeroFace) {
    orientation = 1; orientation[zeroFace] = 0;
    position(zeroFace) = 0.0;
    position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
    position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
    position(zeroFace) = 1.0;
    position(1-zeroFace) = 0.25; spoke.addValue(orientation, position);
    position(1-zeroFace) = 0.75; spoke.addValue(orientation, position);
  }
  Fields_t spokeFlux		(spoke, meshLayout, origin, spacings);

  // Face-centered.

  Centering<Dim> disFace = canonicalCentering<Dim>(FaceType, Discontinuous);


  /* INITIALIZATION */

#ifdef PSEUDOCODE
  // Initialize tensors.
  // Initialize the pressures.
  // Initialize coordinates.
#endif // PSEUDOCODE


  /* COMPUTATION */

#ifdef PSEUDOCODE
  // Compute pressureGradients by simultaneously solving several
  // linear equations.  The operands have different centerings.
  // FIXME
  pressureGradients =
    linearAlgebra<2>(pressure /* cell-centered */,
		     /* Interpolate from vertex-centered to cell-centered: */
		     interpolate<Cell,Vertex>(coordinates),
		     permeability /* cell-centered */,
		     normals /* face-centered */);
#endif // PSEUDOCODE

  // Compute the spoke fluxes.

  // We must multiply three quantities, each with a different
  // centering, to yield values at a fourth-centering.  permeability
  // is cell-centered.  pressureGradient is subcell-centered.  The
  // normals are face-centered.  The product is spoke-centered.

  spokeFlux = 
    dot(replicate(dot(replicate(permeability, nearestNeighbors(cell, subcell)),
		      pressureGradient),
		  nearestNeighbors(subcell, spoke)),
	replicate(meshLayout.unitCoordinateNormals(),
		  nearestNeighbors(disFace, spoke)));

  // Sum the spoke fluxes into a cell flux.

  // Q = \sum_{faces f of cell} sign_f area_f \sum_{subfaces sf of f} q_{sf}
  // We compute this in three steps:
  // 1. Add together the flux values on each face to form a
  //    face-centered field.
  // 2. Multiply each value by the magnitude of the face's normal.
  // 3. Add together each face's value.

  totalFlux =
    sum(spokeFlux.mesh().normals().signedMagnitude() *
	sum(spokeFlux, nearestNeighbors(spoke, disFace)),
	nearestNeighbors(disFace, cell));


  /* TERMINATION */

  std::cout << "total flux:\n" << totalFlux << std::endl;
#ifdef DEBUG
  std::cout << "End program." << std::endl;
#endif // DEBUG
  Pooma::finalize();
  return EXIT_SUCCESS;
}