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