Chevron Code Using New Field Abstractions

[Jeffrey Oldham]

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.
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.


template<class GeometryTag, class T, class Expr, int Dim>
inline
typename Field<GeometryTag, T, Expr>::T_t
faceWeightedSum(const Field<GeometryTag, T, Expr>& inputField, 
		const FieldOffsetList<Dim,Type> &lst,
		const Field<GeometryTag, T, Expr>& outputField)
{
  typedef typename Field<GeometryTag, T, Expr>::T_t T_t;
  typedef typename FieldOffsetList<Dim,Type>::size_type size_type;
  CTAssert((Field<GeometryTag, T, Expr>::dimensions == Dim));

  
	
	      

  // HERE
  
  const size_type lstLength = lst.size();
  PInsist(lstLength > 0, "faceWeightedSum() must be given a nonempty list.");
  T_t init = inputField(lst[0], loc);

  // FIXME inputField.mesh().face(arg).normal() returns a normal
  // vector with length equal to the face's area and direction
  // perpendicular to the face.

  for (size_type i = 1; i < lstLength ; ++i)
    init +=
      outputField.mesh().face(WHICH).normal() *
      inputField(lst[i], loc);
  return init;
}


/** 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)),
	findFieldOffsetList(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;
}