RFC: Hydrodynamics Kernel

Wed Apr 4 18:55:24 UTC 2001

Attached is a partially-completed hydrodynamics kernel program.  Upto
now, I have tried to program clean code, using Pooma abstractions and
trying to determine (and implement) the additional Pooma abstractions
needed to produce clean code.  Given the 15April deadline, soon, I
will begin hacking unclean code to ensure we have working code by the
deadline.

Comments regarding the two conceptual choices at the top of the code
are solicited.

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

// Hydrodynamics Kernel Written Using POOMA

// Following is pseudocode for the hydrodynamics kernel program.  When
// the pseudocode is replaced by Pooma code, the program should
// perform one iteration of a predictor, but not corrector, scheme.
// The program should compile since unimplemented portions are protected by
// "PSEUDOCODE" preprocessor symbols.
// 
// To complete the problem, two conceptual choices remain.
// 1. How should corner masses and corner forces be implemented?
// 
//    Central to the staggered grid concept, a corner is the intersection
//    of two overlapping grids, one shifted half a unit in all dimensions
//    from the other.  For example, in a cell from a two-dimensional (2-D)
//    grid, there are four corner values, which I have marked using
//    2-D binary numbers.
// 
// 			|-------------------|
// 			|01               11|
// 			|         c         |
// 			|00               10|
// 			|-------------------|
// 
//    In three dimensions, there are eight corners in a cell.  Note that
//    referring to a corner value requires both (a cell and a binary
//    number) or (a vertex and a binary number).
// 
//    The pseudocode has ??? operations on these corner fields:
//    1. sumAroundCell(CornerField): For each cell, add together all its
//       corners.  Form a ConstField with the sums at each cell.
//    2. sumAroundVertex(CornerField): For each vertex, add together all
//       its corners.  Form a ConstField with the sums at each vertex.
//    3. computeCornerVolumes(coordinates): Form a CornerField using the
//       coordinates Field as input.  This involves one computation per
//       corner.
//    4. computeCornerForces(pressure, coordinates): Form a CornerField
//       using the two given Fields as input.  This involves one
//       computation per corner.  Pseudocode, which is repetitious, is
//       given.
//    5. Mathematical operation \dot(CornerField, vertex-centered-field).
//       Note the two operands have different number of values, but the
//       "right thing" should happen, i.e., each all corner values around
//       a vertex should each be dotted with the corresponding vertex
//       value.
// 
// 2. Field Stenciling.
// 
//    vertex-centered-field = computeVolumes(cell-centered-field)
// 
//    Scott Haney says that NewField stencils are either broken or
//    difficult (I do not remember which).  Instead he suggested using
//    data-parallel statements, but this requires producing a "parallel"
//    version of the product operation.  Will NewField stencils be fully
//    supported in Pooma 2.4?
// 
// Unfinished Coding Tasks:
// 1. Decide corner field implementation and then implement the
//    operations on it.
// 2. Finish implementing computeVolumes() Field operation.
// 3. Convert to cylindrical or Lagrangian coordinates.
// 4. Improve initialization of "coordinates" Field.
// 5. Improve the implementation of velocity boundary conditions.

// *******************************************************************

#include <iostream>
#include <stdlib.h>
#include <cmath>
#include "Pooma/NewFields.h"

// This hydrodynamics program implements "The Construction of
// Compatible Hydrodynamics Algorithms Utilizing Conservation of Total
// Energy," by E.J. Caramana, D.E. Burton, M.J. Shashkov, and
// P.P. Whalen, \emph{Journal of Computational Physics}, vol. 146
// (1998), pp. 227-262.

// Forward Declarations
template <class Geometry, class VtxT, class T, class Engine>
void computeVolumes(const Field<Geometry,VtxT,Engine> & vtxCentered,
		    Field<Geometry,T,Engine> & output);

template <class Geometry, class VtxT, class Engine>
void enforceZeroPerpendicularComponent(Field<Geometry,VtxT,Engine> & f,
				       const Interval<1> & xInterval,
				       const Interval<1> & yInterval);

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

  // Values
  const double gamma = 4.0/3;		// gas pressure constant $\gamma$
  const double timeStep = 0.01;		// $\Delta t$
  unsigned nuIterations = 1;		// number of iterations

  // Create a simple layout.
  // TODO: Change to Cylindrical coordinates?
  const unsigned Dim = 2;		// Work in a 2D world.
  const unsigned nHorizVerts = 11;	// number of horizontal vertices
  const unsigned nAngles = 5;		// number of angles
  Interval<Dim> vertexDomain;
  vertexDomain[0] = Interval<1>(nHorizVerts);
  vertexDomain[1] = Interval<1>(nAngles);
  DomainLayout<Dim> layout(vertexDomain, GuardLayers<2>(1));

  // Preparation for Field creation.
  Vector<Dim> origin(0.0);
  Vector<Dim> spacings(1.0,1.0);  // TODO: What does this do?
  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 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.
  Cell cell;
  Fields_t internalEnergy  (cell, layout, origin, spacings);
  ConstFields_t zoneMass   (cell, layout, origin, spacings);
  Fields_t pressure	   (cell, layout, origin, spacings);
  Fields_t pressureDensity (cell, layout, origin, spacings);
  Fields_t volume	   (cell, layout, origin, spacings);

  // Vertex-centered Fields.
  Vert vert;
  ConstFields_t pointMass  (vert, layout, origin, spacings);
  Fieldv_t coordinates	   (vert, layout, origin, spacings);
  Fieldv_t velocity	   (vert, layout, origin, spacings);
  Fieldv_t velocityChange  (vert, layout, origin, spacings);

  // Corner Fields.
#ifdef PSEUDOCODE
  // TODO: How should I implement corner Fields?
  Fieldv_t cornerForce	   (ReplicateSubFields<Cell>(1<<Dim), layout,
			    origin, spacings);
  Fields_t cornerMass      (ReplicateSubFields<Cell>(1<<Dim), layout,
			    origin, spacings);
#endif // PSEUDOCODE

  // Initialization
  // TODO: Is coordinates initialization necessary?  What is the best way to do this?
  for (unsigned xIndex = 0; xIndex <= vertexDomain[0].last (); ++xIndex)
    for (unsigned yIndex = 0; yIndex <= vertexDomain[1].last (); ++yIndex)
      coordinates(xIndex,yIndex) = Vector<2>(xIndex,yIndex);
#ifdef DEBUG
  std::cout << "initial coordinates:\n" << coordinates << std::endl;
#endif // DEBUG
  internalEnergy = 1.0;
  pressureDensity = 1.0;
#ifdef PSEUDOCODE
  cornerMass = pressureDensity * computeCornerVolumes(coordinates);
  pointMass = sumAroundCell(cornerMass);
  zoneMass = sumAroundVertex(cornerMass);
#endif // PSEUDOCODE
  velocity = Vector<Dim>(0.0);
  velocityChange.addUpdaters(AllConstantFaceBC<Vector<Dim> >(Vector<Dim>(0.0), false));

  // Iterations
  for (; nuIterations > 0; --nuIterations) {
    pressure = (gamma - 1.0) * pressureDensity * internalEnergy;

#ifdef PSEUDOCODE
    cornerForce = computeCornerForces(pressure, coordinates);
    velocityChange =
      (timeStep / pointMass) * sumAroundCell(cornerForce);
    // TODO: Is there a better way to enforce boundary conditions?
    velocityChange.update();
    enforceZeroPerpendicularComponent(velocityChange,
				      vertexDomain[0], vertexDomain[1]);
    internalEnergy +=
      -(timeStep / pointMass) *
      sumAroundVertex(\dot(cornerForce, velocity + 0.5*velocityChange));
#endif // PSEUDOCODE

    coordinates += (velocity + 0.5*velocityChange) * timeStep;
    velocity += velocityChange;
#ifdef PSEUDOCODE
    volume = computeVolumes(coordinates);
#endif // PSEUDOCODE
    pressureDensity = zoneMass / volume;
  }

  // Termination
  std::cout << "final coordinates:\n" << coordinates << std::endl;

#ifdef DEBUG
  std::cout << "Goodbye, world." << std::endl;
#endif // DEBUG
  Pooma::finalize();
  return EXIT_SUCCESS;
}

// Helper Functions

// This Vector operation proves useful, but I do not know the name of
// the corresponding mathematical operation.
// TODO: Change to specialization for Dim=2.
template <int Dim, class T, class E>
Vector<Dim,T,E>
product(const Vector<Dim,T,E> &v, const Vector<Dim,T,E> &w)
{
  Vector<2,T,E> answer;
  answer(0) = v(0) * w(1);
  answer(1) = -v(1) * w(0);
  return answer;
}

// Compute each cell's volume.
// TODO: Change to handle any dimension or even 2D better.
// TODO: Rewrite using a FieldStencil.
// Geometry: Geometry of the two Field's.
// VtxT: Vector type of Vertex-centered Field
// T: scalar type of Cell-centered Field
// Engine: Engine of the two Field's.
// output should be cell-centered.
template <class Geometry, class VtxT, class T, class Engine>
void computeVolumes(const Field<Geometry,VtxT,Engine> & vtxCentered,
		    Field<Geometry,T,Engine> & output)
{
  // TODO: Check this computation for off-by-one issues.
  Field<Geometry,T,Engine>::Domain_t range = output.physicalCellDomain();
  output(range) = 0.5 *
    std::abs (sum (product (vtxCentered(range+Loc<2>(1,1))-vtxCentered(range+Loc<2>(0,0)),
			    vtxCentered(range+Loc<2>(1,0))-vtxCentered(range+Loc<2>(0,1)))));
  return;
}

// Ensure the field's boundaries' vectors have zero perpendicular
// components.
// TODO: Change to handle any dimension or even 2D better.
template <class Geometry, class VtxT, class Engine>
void enforceZeroPerpendicularComponent(Field<Geometry,VtxT,Engine> & f,
				       const Interval<1> & xInterval,
				       const Interval<1> & yInterval)
{
  f(0,yInterval).comp(0) = 0.0;
  f(xInterval,0).comp(1) = 0.0;
  return;
}

#ifdef PSEUDOCODE
// Compute the corner forces.
computeCornerForces(const Field & pressure, const Field & coordinates,
		    Field & cornerForces)
{
  Field<Geometry2,T,Engine>::Domain_t range = input.physicalCellDomain();
  Loc<Dim> previous_loc, loc, next_loc;
  previous_loc = Loc<2>(1,0); loc = Loc<2>(1,1); next_loc = Loc<2>(0,1);
  cornerForces[zoneToBinaryCorner(loc)](range) =
    pressure(range) * 0.5 *
    product(Vector<Dim>(1.0),
	    coordinates(range+next_loc) - coordinates(range+previous_loc));
  previous_loc = loc; loc = next_loc; next_loc = Loc<2>(0,0);
  cornerForces[zoneToBinaryCorner(loc)](range) =
    pressure(range) * 0.5 *
    product(Vector<Dim>(1.0),
	    coordinates(range+next_loc) - coordinates(range+previous_loc));
  previous_loc = loc; loc = next_loc; next_loc = Loc<2>(1,0);
  cornerForces[zoneToBinaryCorner(loc)](range) =
    pressure(range) * 0.5 *
    product(Vector<Dim>(1.0),
	    coordinates(range+next_loc) - coordinates(range+previous_loc));
  previous_loc = loc; loc = next_loc; next_loc = Loc<2>(1,1);
  cornerForces[zoneToBinaryCorner(loc)](range) =
    pressure(range) * 0.5 *
    product(Vector<Dim>(1.0),
	    coordinates(range+next_loc) - coordinates(range+previous_loc));
  return;
}
#endif // PSEUDOCODE