EulerUtil.H

Go to the documentation of this file.

#ifndef __EULER_UTIL_H__
#define __EULER_UTIL_H__

#include "PlasCom2.H"
#include "PCPPCommUtil.H"
#include "PCPPIntervalUtils.H"
#include "PCPPReport.H"

namespace euler {
  namespace util {

/*
    // MJA, to be deleted
    // Input:
    // inInterval:  (xStart xEnd yStart yEnd zStart zEnd)
    // bufferSize:  (numX numY numZ)
    // stateBuffer: [ *rho *rhoV *rhoE *scalarVars ]
    //
    // Output:
    // dvBuffer:    [1/rho p T v1 v2 v3 ]
    int ComputeDVBufferIdeal(const pcpp::IndexIntervalType &regionInterval,
                        const std::vector<size_t>     &bufferSizes,
                        const std::vector<double *>   &stateBuffers,
                        std::vector<double *>         &dvBuffers,
                        const double                   gamma,
                        const double                   specGasConst);
*/

    // Input:
    // inInterval:  (xStart xEnd yStart yEnd zStart zEnd)
    // bufferSize:  (numX numY numZ)
    // stateBuffer: [ *rho *rhoV *rhoE *scalarVars ]
    //
    // Output:
    // dvBuffer:    [p T 1/rho v1 v2 v3 rho*e]
    int ComputeDVBuffer(const pcpp::IndexIntervalType &regionInterval,
                        const std::vector<size_t>     &bufferSizes,
                        const std::vector<double *>   &stateBuffers,
                        std::vector<double *>         &dvBuffers);
    
    

    // Fortran-like interface 
    int ComputeDVBuffer2(const int *numDimPtr,
                        const size_t *bufferSize,
                        const size_t *numPointsPtr,
                        const size_t *bufferInterval,
                        const double *rhoBuffer,
                        const double *rhoVBuffer,
                        const double *rhoEBuffer,
                        double *pressureBuffer,
                        double *tempBuffer,
                        double *rhom1Buffer,
                        double *velBuffer);

    inline double Pulse(double amp,double width,double centerX,
                        double centerY,double centerZ,double x,double y,double z)
    {
      double xTerm     = (x-centerX);
      double yTerm     = (y-centerY);
      double zTerm     = (z-centerZ);
      double eTerm     = zTerm*zTerm + xTerm*xTerm + yTerm*yTerm;
      return(amp*std::exp(-(eTerm)/width));  
    };

    inline double Window(double tWidth,double wWidth,
                         double centerX,double centerY,double centerZ,
                         double x,double y,double z)
    {

      double xWindow   = 1.0;
      double xPos = std::sqrt((x-centerX)*(x-centerX)+(y-centerY)*(y-centerY)+
                              (z-centerZ)*(z-centerZ));
      double A = 3.14159265358979*2.0/tWidth;
      double wPos = wWidth - xPos;
      xWindow *= (std::tanh(A*wPos)+1.0)/2.0;
      return(xWindow);
    };

    inline double Parabola(double amp,double centerX,
                           double centerY,double centerZ,double x,double y,double z)
    {
      double xTerm     = (x-centerX);
      double yTerm     = (y-centerY);
      double zTerm     = (z-centerZ);
      double eTerm     = zTerm*zTerm + xTerm*xTerm + yTerm*yTerm;
      return(amp*eTerm);  
    };

    template<typename GridType,typename HaloType>
    int CreateSimulationFixtures(GridType &inGrid,HaloType &inHalo,
                                 plascom2::operators::sbp::base &inOperator,
                                 const std::vector<size_t> &gridSizes,
                                 int interiorOrder,
                                 pcpp::CommunicatorType &inComm,
                                 std::ostream *messageStreamPtr = NULL)
    {
  
      int numProc = inComm.Size();
      int myRank  = inComm.Rank();
  
      if(!messageStreamPtr)
        messageStreamPtr = &std::cout;

      pcpp::CommunicatorType &gridComm(inGrid.Communicator());

      int numDim = gridSizes.size();
      std::vector<double> dX(numDim,0);
      size_t numGlobalNodes = 1;
      size_t numGlobalCells = 1;
  
      for(int iDim = 0;iDim < numDim;iDim++){
        numGlobalNodes *= gridSizes[iDim];
        numGlobalCells *= (gridSizes[iDim]-1);
        dX[iDim] = 1.0/static_cast<double>(gridSizes[iDim]-1);
      }

      inGrid.SetGridSizes(gridSizes);
      inGrid.SetGridSpacings(dX);
      
      plascom2::operators::sbp::Initialize(inOperator,interiorOrder);
      //      int boundaryDepth = (inOperator.numStencils-1)/2;
      int boundaryDepth = inOperator.boundaryDepth;
      std::vector<int> haloDepths(2*numDim,boundaryDepth);

      pcpp::ParallelTopologyInfoType cartInfo;
      cartInfo.numDimensions = numDim;
      pcpp::ParallelTopologyInfoType pbsCartInfo;
      pbsCartInfo.numDimensions = numDim;
      pbsCartInfo.isPeriodic.resize(numDim,1);
      pcpp::comm::SetupCartesianTopology(inComm,pbsCartInfo,gridComm,*messageStreamPtr);

      // Grab the local coordinates and global dimension of the Cartesian topo
      std::vector<int> &cartCoords(gridComm.CartCoordinates());
      std::vector<int> &cartDims(gridComm.CartDimensions());

      // Generate a report of the Cartesian setup and put it the messageStream
      pcpp::report::CartesianSetup(*messageStreamPtr,pbsCartInfo);

      pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      pcpp::IndexIntervalType globalInterval;
      globalInterval.InitSimple(gridSizes);

      // === Partition the grid ===
      if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,
                                                partitionInterval,*messageStreamPtr)){
        *messageStreamPtr << "TestFixtures:SetupSimulationFixtures:Error: Partitioning failed." << std::endl;
        return(1);
      }

      std::vector<size_t> extendGrid(2*numDim,0);
      std::vector<bool> isPeriodic(numDim,true);
      size_t numPointsPart = partitionInterval.NNodes();
  
      for(int iDim = 0;iDim < numDim; iDim++){
        if(isPeriodic[iDim]){
          extendGrid[2*iDim]   = haloDepths[2*iDim];
          extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
        } else {
          if(partitionInterval[iDim].first == 0)
            extendGrid[2*iDim] = haloDepths[2*iDim];
          if(partitionInterval[iDim].second == (gridSizes[iDim]-1))
            extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
        }
      }

      inGrid.SetDimensionExtensions(extendGrid);
      if(inGrid.Finalize()){
        *messageStreamPtr << "TestFixtures:SetupSimulationFixtures:Error: Grid finalization failed." << std::endl;
        return(1);
      }

      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());  
      size_t numPointsBuffer = inGrid.BufferSize();

      // Set up the HALO data structure 
      //  halo_t testHalo(globalInterval,pbsPartInterval);
      inHalo.SetGridInterval(globalInterval);
      inHalo.SetPartitionInterval(partitionInterval);

      // This forces the halo to have neighbors even on 
      // partitions with domain boundaries.
      std::vector<bool> haveNeighbors(2*numDim,true);
      inHalo.SetNeighbors(haveNeighbors);

      std::vector<pcpp::IndexIntervalType> remoteHaloExtents(inHalo.CreateRemoteHaloExtents(extendGrid));
      std::vector<pcpp::IndexIntervalType> localHaloExtents(inHalo.CreateLocalHaloExtents(extendGrid));
  
      inHalo.SetLocalBufferSizes(bufferSizes);
      inHalo.SetLocalPartitionExtent(partitionBufferInterval);
  
      inHalo.SetRemoteHaloExtents(remoteHaloExtents);
      inHalo.SetLocalHaloExtents(localHaloExtents);
  
      const std::vector<pcpp::IndexIntervalType> &remoteHaloBufferExtents(inHalo.RemoteHaloBufferExtents());
      const std::vector<pcpp::IndexIntervalType> &localHaloBufferExtents(inHalo.LocalHaloBufferExtents());
  
      inHalo.CreateSimpleSendIndices();
      inHalo.CreateSimpleRecvIndices();
  
      return(0);

    };
 
    //     template<typename GridType>
    //     int CreateSimulationFixtures(GridType &inGrid,
    //                                  plascom2::operators::sbp::base &inOperator,
    //                                  std::vector<double> &physicalExtent,
    //                                  const std::vector<size_t> &gridSizes,
    //                                  int interiorOrder,
    //                                  pcpp::CommunicatorType &inComm,
    //                                  std::ostream *messageStreamPtr = NULL)

    template<typename GridType>
    int InitializeSimulationFixtures(GridType &inGrid,
                                     plascom2::operators::sbp::base &inOperator,
                                     int interiorOrder,
                                     pcpp::CommunicatorType &inComm,
                                     std::ostream *messageStreamPtr = NULL)
    {
      typedef typename GridType::HaloType HaloType;

      int numProc = inComm.Size();
      int myRank  = inComm.Rank();
  
      if(!messageStreamPtr)
        messageStreamPtr = &std::cout;

      pcpp::CommunicatorType &gridComm(inGrid.Communicator());

      std::vector<size_t> &gridSizes(inGrid.GridSizes());
      std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
      std::vector<bool>   &isPeriodic(inGrid.PeriodicDirs());

      int numDim = gridSizes.size();
      std::vector<double> dX(numDim,0);
      size_t numGlobalNodes = 1;
      size_t numGlobalCells = 1;
  
      if(physicalExtent.empty()){
        physicalExtent.resize(numDim*2,0);
        for(int iDim = 0;iDim < numDim;iDim++)
          physicalExtent[2*iDim+1] = 1.0;
      }

      if(isPeriodic.empty()){
        isPeriodic.resize(numDim,true);
      }

      for(int iDim = 0;iDim < numDim;iDim++){
        numGlobalNodes *= gridSizes[iDim];
        numGlobalCells *= (gridSizes[iDim]-1);
        dX[iDim] = (physicalExtent[2*iDim+1]-physicalExtent[2*iDim])/static_cast<double>(gridSizes[iDim]-1);
      }

      inGrid.SetGridSpacings(dX);
      
      plascom2::operators::sbp::Initialize(inOperator,interiorOrder);

      int boundaryDepth = inOperator.boundaryDepth;
      std::vector<int> haloDepths(2*numDim,boundaryDepth);

      pcpp::ParallelTopologyInfoType cartInfo;
      cartInfo.numDimensions = numDim;
      pcpp::ParallelTopologyInfoType pbsCartInfo;
      pbsCartInfo.numDimensions = numDim;
      pbsCartInfo.isPeriodic.resize(numDim,0);
      for(int iDim = 0;iDim < numDim;iDim++)
        pbsCartInfo.isPeriodic[iDim] = isPeriodic[iDim];
      
      pcpp::comm::SetupCartesianTopology(inComm,pbsCartInfo,gridComm,*messageStreamPtr);

      // Grab the local coordinates and global dimension of the Cartesian topo
      std::vector<int> &cartCoords(gridComm.CartCoordinates());
      std::vector<int> &cartDims(gridComm.CartDimensions());

      std::vector<int> gridNeighbors;
      gridComm.CartNeighbors(gridNeighbors);


      // Generate a report of the Cartesian setup and put it the messageStream
      pcpp::report::CartesianSetup(*messageStreamPtr,pbsCartInfo);

      pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      pcpp::IndexIntervalType globalInterval;
      globalInterval.InitSimple(gridSizes);

      // === Partition the grid ===
      if(pcpp::util::PartitionCartesianInterval(globalInterval,cartDims,cartCoords,
                                                partitionInterval,*messageStreamPtr)){
        *messageStreamPtr << "TestFixtures:SetupSimulationFixtures:Error: Partitioning failed." << std::endl;
        return(1);
      }

      std::vector<size_t> extendGrid(2*numDim,0);
      //      std::vector<bool> isPeriodic(numDim,true);
      size_t numPointsPart = partitionInterval.NNodes();
  
      std::vector<bool> haveNeighbors(2*numDim,true);
      for(int iDim = 0;iDim < numDim; iDim++){
        if(isPeriodic[iDim]){
          extendGrid[2*iDim]   = haloDepths[2*iDim];
          extendGrid[2*iDim+1] = haloDepths[2*iDim+1];
        } else {
          if(partitionInterval[iDim].first == 0) {
            haveNeighbors[2*iDim] = false;
          } else { 
            extendGrid[2*iDim] = haloDepths[2*iDim];
          }
          if(partitionInterval[iDim].second == (gridSizes[iDim]-1)){
            haveNeighbors[2*iDim+1] = false;
          } else {
            extendGrid[2*iDim + 1] = haloDepths[2*iDim+1];
          }
        }
      }

      inGrid.SetDimensionExtensions(extendGrid);


      if(inGrid.Finalize()){
        *messageStreamPtr << "TestFixtures:SetupSimulationFixtures:Error: Grid finalization failed." << std::endl;
        return(1);
      }
      
      if(inGrid.GenerateCoordinates(*messageStreamPtr)){
        *messageStreamPtr << "TestFixtures:SetupSimulationFixtures:Error: Grid coordinates generation failed." 
                          << std::endl;
        return(1);
      }
      
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());  
      size_t numPointsBuffer = inGrid.BufferSize();

      // Set up the HALO data structure 
      //  halo_t testHalo(globalInterval,pbsPartInterval);
      HaloType &inHalo(inGrid.Halo());
      inHalo.SetGridInterval(globalInterval);
      inHalo.SetPartitionInterval(partitionInterval);
      inHalo.SetPeriodicDirs(pbsCartInfo.isPeriodic);

      // This forces the halo to have neighbors even on 
      // partitions with domain boundaries.
      int rank = gridComm.Rank();
      int numProcGrid = gridComm.NProc();

      inHalo.SetNeighbors(haveNeighbors);

      std::vector<pcpp::IndexIntervalType> remoteHaloExtents(inHalo.CreateRemoteHaloExtents(extendGrid));
      std::vector<pcpp::IndexIntervalType> localHaloExtents(inHalo.CreateLocalHaloExtents(extendGrid));

      inHalo.SetLocalBufferSizes(bufferSizes);
      inHalo.SetLocalPartitionExtent(partitionBufferInterval);
  
      inHalo.SetRemoteHaloExtents(remoteHaloExtents);
      inHalo.SetLocalHaloExtents(localHaloExtents);
  
      const std::vector<pcpp::IndexIntervalType> &remoteHaloBufferExtents(inHalo.RemoteHaloBufferExtents());
      const std::vector<pcpp::IndexIntervalType> &localHaloBufferExtents(inHalo.LocalHaloBufferExtents());
  
      inHalo.CreateSimpleSendIndices();
      inHalo.CreateSimpleRecvIndices();
  
      return(0);

    };

    template<typename GridType,typename StateType>
    int SetupEulerState(const GridType &inGrid,StateType &inState,StateType &inParams,int numScalars,
                        bool withFlux = false)
    {

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<double> &dX(inGrid.DX());
      int numDim = bufferSizes.size();

      inState.AddField("simTime",'s',1,8,"s");
      
      // Conserved quantities, mass, momentum, and energy densities
      inState.AddField("rho",'n',1,8,"mass/vol");
      inState.AddField("rhoV",'n',numDim,8,"momentum/vol");
      inState.AddField("rhoE",'n',1,8,"energy/vol");

      if(numScalars>0)
        inState.AddField("scalarVars",'n',numScalars,8,"m/M");
      
      //for(int iScalar = 0;iScalar < numScalar;iScalar++){
        //std::ostringstream Ostr;
        //Ostr << "scalar-" << iScalar+1;
        //inState.AddField(Ostr.str(),'n',1,8,"m/M");
      //}

      // Dependent quantities, pressure, temperature, velocity
      inState.AddField("pressure",'n',1,8,"pressure");
      inState.AddField("temperature",'n',1,8,"temperature");
      inState.AddField("rhom1",'n',1,8,"volume");
      inState.AddField("velocity",'n',numDim,8,"velocity");

      // Flux quanitites (if tracking desired)
      if(withFlux){
        inState.AddField("massFlux",'n',numDim,8,"");
        inState.AddField("mom1Flux",'n',numDim,8,"");
        inState.AddField("mom2Flux",'n',numDim,8,"");
        if(numDim == 3){
          inState.AddField("mom3Flux",'n',3,8,"");
        }
        inState.AddField("energyFlux",'n',3,8,"");
        if(numScalars>0)
          inState.AddField("scalarVarsFlux",'n',numScalars,8,"");
      }

      inState.Create(numPointsBuffer,0);
      if(numScalars == 0)
        inState.SetStateFields("rho rhoV rhoE");
      else
        inState.SetStateFields("rho rhoV rhoE scalarVars");
        
      
      inParams.AddField("gamma",'s',1,8,"");
      inParams.AddField("inputDT",'s',1,8,"s");
      inParams.AddField("inputCFL",'s',1,8,"");
      //inParams.AddField("refC",'s',1,8,"soundspeed");
      //inParams.AddField("refR",'s',1,8,"gasconstant");
      inParams.AddField("refRe",'s',1,8,"ReynoldsNo");
      inParams.AddField("refPr",'s',1,8,"PrandtlNo");

      inParams.Create(numPointsBuffer,0);
     
      return(0);
    };

    template<typename DatasetType>
    int ValidateState(const DatasetType &inState,
                      const DatasetType &inParam,
                      std::ostream &outStream)
    {
      int RHO     = 1<<1;
      int RHOV    = 1<<2;
      int RHOE    = 1<<3;
      int SIMT    = 1<<4;
      int INDT    = 1<<5;
      int INCFL   = 1<<6;
      int REFL    = 1<<7;
      int REFP    = 1<<8;
      int REFT    = 1<<9;
      int REFRE   = 1<<13;
      int REFRHO  = 1<<10;
      int GAMMA   = 1<<11;
      int NOND    = 1<<12;

      int dataFaults = 0;

      int dataHandle = inState.GetDataIndex("rho");
      if(dataHandle < 0){
        dataFaults |= RHO;
      }
      dataHandle = inState.GetDataIndex("rhoV");
      if(dataHandle < 0){
        dataFaults |= RHOV;
      }
      dataHandle = inState.GetDataIndex("rhoE");
      if(dataHandle < 0){
        dataFaults |= RHOE;
      }
      dataHandle = inState.GetDataIndex("simTime");
      if(dataHandle < 0){
        dataFaults |= SIMT;
      }
      dataHandle = inParam.GetDataIndex("inputDT");
      if(dataHandle < 0){
        dataFaults |= INDT;
      }
      dataHandle = inParam.GetDataIndex("inputCFL");
      if(dataHandle < 0){
        dataFaults |= INCFL;
      }
      dataHandle = inParam.GetDataIndex("refRho");
      if(dataHandle < 0){
        dataFaults |= REFRHO;
      }
      dataHandle = inParam.GetDataIndex("refPressure");
      if(dataHandle < 0){
        dataFaults |= REFP;
      } 
      dataHandle = inParam.GetDataIndex("refTemperature");
      if(dataHandle < 0){
        dataFaults |= REFT;
      }
      dataHandle = inParam.GetDataIndex("refRe");
      if(dataHandle < 0){
        dataFaults |= REFRE;
      }
      dataHandle = inParam.GetDataIndex("gamma");
      if(dataHandle < 0){
        dataFaults |= GAMMA;
      }
      dataHandle = inParam.GetDataIndex("refLength");
      if(dataHandle < 0){
        dataFaults |= REFL;
      }
      dataHandle = inParam.GetDataIndex("nonDimensional");
      if(dataHandle < 0){
        dataFaults |= NOND;
      }
     
      int errorState = 0;
      if(dataFaults&RHO){
        outStream << "Required data \"rho\" did not exist in state." << std::endl;
        errorState++;
      }
      if(dataFaults&RHOV){
        outStream << "Required data \"rhoV\" did not exist in state." << std::endl;
        errorState++;
      }
      if(dataFaults&RHOE){
        outStream << "Required data \"rhoE\" did not exist in state." << std::endl;
        errorState++;
      }
      if(dataFaults&SIMT){
        outStream << "Required data \"simTime\" did not exist in state." << std::endl;
        errorState++;
      }
      if(dataFaults&INDT){
        outStream << "Parameter data \"inputDT\" did not exist." << std::endl;
      }
      if(dataFaults&INCFL){
        outStream << "Parameter data \"inputCFL\" did not exist." << std::endl;
      }
      if(dataFaults&REFL){
        outStream << "Parameter data \"refLength\" did not exist." << std::endl;
      }
      if(dataFaults&REFP){
        outStream << "Parameter data \"refPressure\" did not exist." << std::endl;
      }
      if(dataFaults&REFT){
        outStream << "Parameter data \"refTemperature\" did not exist." << std::endl;
      }
      if(dataFaults&REFRHO){
        outStream << "Parameter data \"refRho\" did not exist." << std::endl;
      }
      if(dataFaults&GAMMA){
        outStream << "Parameter data \"gamma\" did not exist." << std::endl;
      }
      if(dataFaults&NOND){
        outStream << "Parameter data \"nonDimensional\" did not exist." << std::endl;
      }
      if(dataFaults&REFRE){
        outStream << "Parameter data \"refRe\" did not exist." << std::endl;
      }
      return(errorState);
    };
    
    template<typename GridType,typename StateType>
    int InitializeQuiescentState(const GridType &inGrid,StateType &inState,
                                 StateType &inParam,int threadId,std::ostream *messageStream = NULL)
    {
      if(!messageStream)
        messageStream = &std::cout;
      
      std::ostringstream Ostr;
      if(ValidateState(inState,inParam,Ostr)){
        std::cout << "Input State in Error:" << std::endl;
        *messageStream << Ostr.str() << std::endl;
        return(1);
      } else {
        *messageStream << "Input State Validation: " << std::endl
                       << Ostr.str() << std::endl;
      }
      
      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      
      int numDim = gridSizes.size();
      
      double *timePtr       = NULL;
      double *deltaTPtr     = NULL;
      double *gammaPtr      = NULL;
      double *CFLPtr        = NULL;
      double *lengthRefPtr  = NULL;
      double *rhoRefPtr     = NULL;
      double *presRefPtr    = NULL;
      double *tempRefPtr    = NULL;
      double *sigmaPtr      = NULL;
      double *rePtr         = NULL;
      int *nonDimenPtr      = NULL;

      int dataHandle = inState.GetDataIndex("simTime");
      if(dataHandle < 0){
        inState.AddField("simTime",'d',1,8,"time");
        inState.Create();
      }
      timePtr       = inState.template GetFieldBuffer<double>("simTime");
      if(*timePtr < 0)
        *timePtr = 0.0;
      
      dataHandle = inParam.GetDataIndex("inputDT");
      if(dataHandle < 0){
        inParam.AddField("inputDT",'d',1,8,"time");
        inParam.Create();
      }
      deltaTPtr     = inParam.template GetFieldBuffer<double>("inputDT");

      dataHandle = inParam.GetDataIndex("gamma");
      if(dataHandle < 0){
        inParam.AddField("gamma",'d',1,8,"");
        inParam.Create();
      }
      gammaPtr     = inParam.template GetFieldBuffer<double>("gamma");

      dataHandle = inParam.GetDataIndex("inputCFL");
      if(dataHandle < 0){
        inParam.AddField("inputCFL",'d',1,8,"");
        inParam.Create();
      }
      CFLPtr     = inParam.template GetFieldBuffer<double>("inputCFL");

      dataHandle = inParam.GetDataIndex("refLength");
      if(dataHandle < 0){
        inParam.AddField("refLength",'d',1,8,"length");
        inParam.Create();
      }
      lengthRefPtr     = inParam.template GetFieldBuffer<double>("refLength");

      dataHandle = inParam.GetDataIndex("refRho");
      if(dataHandle < 0){
        inParam.AddField("refRho",'d',1,8,"density");
        inParam.Create();
      }
      rhoRefPtr     = inParam.template GetFieldBuffer<double>("refRho");

      dataHandle = inParam.GetDataIndex("refPressure");
      if(dataHandle < 0){
        inParam.AddField("refPressure",'d',1,8,"pressure");
        inParam.Create();
      }
      presRefPtr     = inParam.template GetFieldBuffer<double>("refPressure");

      dataHandle = inParam.GetDataIndex("refTemperature");
      if(dataHandle < 0){
        inParam.AddField("refTemperature",'d',1,8,"temperature");
        inParam.Create();
      }
      tempRefPtr     = inParam.template GetFieldBuffer<double>("refTemperature");

      dataHandle = inParam.GetDataIndex("sigma");
      if(dataHandle < 0){
        inParam.AddField("sigma",'d',1,8,"");
        inParam.Create();
      }
      sigmaPtr     = inParam.template GetFieldBuffer<double>("sigma");

      dataHandle = inParam.GetDataIndex("refRe");
      if(dataHandle < 0){
        inParam.AddField("refRe",'d',1,8,"ReynoldsNo");
        inParam.Create();
      }
      rePtr     = inParam.template GetFieldBuffer<double>("refRe");

      dataHandle = inParam.GetDataIndex("nonDimensional");
      if(dataHandle < 0){
        inParam.AddField("nonDimensional",'d',1,4,"");
        inParam.Create();
      }
      nonDimenPtr    = inParam.template GetFieldBuffer<int>("nonDimensional");

      *messageStream << "Parameters:    "                    << std::endl
                     << "time:           " << *timePtr        << std::endl
                     << "deltaT:         " << *deltaTPtr      << std::endl
                     << "gamma:          " << *gammaPtr       << std::endl
                     << "CFL:            " << *CFLPtr         << std::endl
                     << "lengthRef:      " << *lengthRefPtr   << std::endl
                     << "rhoRef:         " << *rhoRefPtr      << std::endl
                     << "pressureRef:    " << *presRefPtr     << std::endl
                     << "temperatureRef: " << *tempRefPtr     << std::endl
                     << "Re:             " << *rePtr          << std::endl
                     << "Sigma:          " << *sigmaPtr       << std::endl
                     << "nonDimensional: " << *nonDimenPtr    << std::endl;
      
      double gammaValue = 1.4;

      // check the reference state. Default to a dimensional mode if set incorrectly.
      if((*presRefPtr < 0 || *presRefPtr > 1e10) ||
         (*rhoRefPtr < 0  || *rhoRefPtr > 1e10)  ||
         (*tempRefPtr < 0 || *tempRefPtr > 1e10)){
        *messageStream << "Invalid reference state in Initialize quiescent state."
                       << " Defaulting to non-dimensional mode." << std::endl;
        *presRefPtr = 101325.0;
        *rhoRefPtr = 1.225;
        *tempRefPtr = 288.15;

        *messageStream << "Reference state: "                << std::endl
                       << "\t Density = "     << *rhoRefPtr  << std::endl
                       << "\t Pressure = "    << *presRefPtr << std::endl
                       << "\t Temperature = " << *tempRefPtr << std::endl;
      }
      
      if(*CFLPtr < 0 || *CFLPtr > 4)
        *CFLPtr    = 1.0;
      
      *timePtr   = 0.0;
      
      if(*gammaPtr < 0 || *gammaPtr > 4)
        *gammaPtr  = gammaValue;
      
      gammaValue = *gammaPtr;
      double cRef = sqrt(gammaValue*(*presRefPtr)/(*rhoRefPtr));

      if(*deltaTPtr < 1e-20 || *deltaTPtr > 2){
        *messageStream << "DT parameter outside bounds (" << *deltaTPtr 
                       << "). Calculating max timestep based on CFL ("
                       << *CFLPtr << ") and reference speed (" << cRef
                       << ")" << std::endl;
        // Set params in here for now
        double minSpacing = 100000.0;
        std::vector<double>::const_iterator dxIt = dX.begin();
        while(dxIt != dX.end()){
          if(*dxIt < minSpacing) minSpacing = *dxIt;
          dxIt++;
        }
        *deltaTPtr = minSpacing*(*CFLPtr)/(cRef);
      }

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      std::vector<double> rhoV(numDim,0.0); // really just velocity at this point
      // Comment next line for quiescent,  uncomment for constant flow in X
      //      rhoV[0] = 1.0;
      double rho       = 1.0;
      double pressure  = 1.0;
      if(!*nonDimenPtr){
        rho       = *rhoRefPtr;
        pressure  = *presRefPtr;
      }
      double eInternal = (pressure/(gammaValue-1.0))/rho;
      double eKinetic  = 0.0;
      std::vector<double>::iterator rhoVIt = rhoV.begin();
      while(rhoVIt != rhoV.end()){
        eKinetic += (rho * (*rhoVIt) * (*rhoVIt));
        *rhoVIt *= rho; // make it rho*V for real
        rhoVIt++;
      }
      eKinetic *= .5;
      double rhoE  = rho*(eInternal + eKinetic);
      
      if(numDim == 3){
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        size_t kStart = partitionBufferInterval[2].first;
        size_t kEnd   = partitionBufferInterval[2].second;
        size_t nPlane = bufferSizes[0]*bufferSizes[1];
        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              rhoPtr[bufferIndex]  = rho;
              rhoEPtr[bufferIndex] = rhoE;
              rhoVPtr[bufferIndex] = rhoV[0];
              rhoVPtr[bufferIndex+numPointsBuffer] = rhoV[1];
              rhoVPtr[bufferIndex+2*numPointsBuffer] = rhoV[2];
            }
          }
        }
      } else if (numDim == 2) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = rhoE;
            rhoVPtr[bufferIndex] = rhoV[0];
            rhoVPtr[bufferIndex+numPointsBuffer] = rhoV[1];
          }
        }
      } else if (numDim == 1) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          rhoPtr[iIndex]  = rho;
          rhoEPtr[iIndex] = rhoE;
          rhoVPtr[iIndex] = rhoV[0];
        }
      }
      
      return(0);
    }

    template<typename GridType,typename StateType>
    int InitializeShock1D(const GridType &inGrid,StateType &inState,
                          StateType &inParamState,const std::vector<double> &inParams, int threadId,
                          std::ostream *messageStream)
    {
      if(!messageStream)
        messageStream = &std::cout;

      double inletX = 0.0; // inlet x coordinate
      double inletY = 0.0; // inlet y coordinate
      double inletZ = 0.0; // inlet z coordinate
      int direction = 1;
      double location = 0.1;
      double transitionWidth = 0.0;
      double mach = 1.0; // shock pressure ratio
      int referenceFrame = 0; // 0 for stationary (lab frame), 1 for moving with shock

      int numParams = inParams.size();
      if(numParams > 0)
        mach = inParams[0];
      if(numParams > 1)
        direction = static_cast<int>(inParams[1]);
      if(numParams > 2)
        location = inParams[2];
      if(numParams > 3)
        transitionWidth = inParams[3];
      if(numParams > 4)
        referenceFrame = inParams[4];
     
      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<double> &coordinateData(inGrid.CoordinateData());
      
      int numDim = gridSizes.size();
      
      double *gammaPtr      = inParamState.template GetFieldBuffer<double>("gamma");
      double *lengthRefPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *rhoRefPtr     = inParamState.template GetFieldBuffer<double>("refRho");
      double *presRefPtr    = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr    = inParamState.template GetFieldBuffer<double>("refTemperature");
      double *rePtr         = inParamState.template GetFieldBuffer<double>("refRe");
      double *prPtr         = inParamState.template GetFieldBuffer<double>("refPr");
      int *nonDimenPtr      = inParamState.template GetFieldBuffer<int>("nonDimensional");

      double gamma = *gammaPtr;
      double refLength = *lengthRefPtr;
      double refRho = *rhoRefPtr;
      double refPressure = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef = sqrt(gamma*refPressure/refRho);
      double velRef = sqrt(refPressure/refRho);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }

      double xMin = readGridExtent[0];
      double xMax = readGridExtent[1];
      double yMin = readGridExtent[2];
      double yMax = readGridExtent[3];
      double zMin = 0.;
      double zMax = 0.;
      if(numDim > 2 ) {
        zMin = readGridExtent[4];
        zMax = readGridExtent[5];
      }


      // dimensional flow conditions
      // calculated from normal shock relations
      //
      double pressureRatio = (2.*gamma*mach*mach-(gamma-1.))/(gamma+1.);
      double densityRatio = (gamma+1.)*mach*mach/((gamma-1.)*mach*mach+2.);
      double mach2 = sqrt(((gamma-1.)*mach*mach+2.)/(2.*gamma*mach*mach-(gamma-1.)));
      double rho1       = *rhoRefPtr;
      double pressure1  = *presRefPtr;
      double rho2 = rho1*densityRatio;
      double pressure2 = pressure1*pressureRatio;
      // convective velocity (behind shock wave in lab frame, shock moving)
      double velocity1 = 0.;
      double velocity2 = -mach*sndSpdRef*(1/densityRatio-1);
      // moving shock reference frame
      if(referenceFrame == 1) {
        velocity1 = -mach*sndSpdRef;
        velocity2 = velocity1/densityRatio;
      }
      
      std::vector<double> rhoV(numDim,0.0);

      // scale for nonDimensionalization
      if(*nonDimenPtr == 1){
        rho1       = rho1/(*rhoRefPtr);
        rho2       = rho2/(*rhoRefPtr);
        pressure1  = pressure1/(*presRefPtr);
        pressure2  = pressure2/(*presRefPtr);
        velocity1 = velocity1/velRef;
        velocity2 = velocity2/velRef;
      } else if(*nonDimenPtr == 2){
        rho1       = rho1/(*rhoRefPtr);
        rho2       = rho2/(*rhoRefPtr);
        pressure1  = *presRefPtr/((*rhoRefPtr)*sndSpdRef*sndSpdRef);
        pressure2  = pressure2/((*rhoRefPtr)*sndSpdRef*sndSpdRef);
        velocity1 = velocity1/sndSpdRef;
        velocity2 = velocity2/sndSpdRef;
      }
      // internal energy
      double rhoE1 = pressure1/(gamma-1.0) + 0.5*rho1*velocity1*velocity1;
      double rhoE2 = pressure2/(gamma-1.0) + 0.5*rho2*velocity2*velocity2;

      if(referenceFrame == 0){
        *messageStream << "Shock moving (lab) reference frame" << std::endl;
      } else if (referenceFrame == 1) {
        *messageStream << "Shock stationary (shock attached) reference frame" << std::endl;
      } else  {
        *messageStream << "Invalid reference frame selection. " 
                  << "Must be 0 (lab reference) or 1 (shock attached)." << std::endl;
        return(1);
      }


      *messageStream << "InitializeShock1D: Pressure Ratio " << pressureRatio << std::endl;
      *messageStream << "InitializeShock1D: Density Ratio " << densityRatio << std::endl;
      *messageStream << "InitializeShock1D: Upstream Mach Number " << mach << std::endl;
      *messageStream << "InitializeShock1D: Upstream Total Energy " 
                << rhoE1+0.5*rho1*densityRatio*velocity1*velocity1 << std::endl;
      *messageStream << "InitializeShock1D: Post-shock Density " << rho2 << std::endl;
      *messageStream << "InitializeShock1D: Post-shock Pressure " << pressure2 << std::endl;
      *messageStream << "InitializeShock1D: Post-shock Velocity " << velocity2 << std::endl;
      *messageStream << "InitializeShock1D: Post-shock Total Energy " 
                << rhoE1*pressureRatio+0.5*rho2*velocity2*velocity2 << std::endl;

      double transitionPoint = 0.2;
      if(direction == 1) {
        transitionPoint = location;
      } else if(direction == 2) {
        transitionPoint = location;
      } else if(direction == 3) {
        transitionPoint = location;
      }

      if(numDim == 3){
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        size_t kStart = partitionBufferInterval[2].first;
        size_t kEnd   = partitionBufferInterval[2].second;
        size_t nPlane = bufferSizes[0]*bufferSizes[1];
        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          size_t kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          double z = kPartIndex*dX[2];
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
            size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
            double y = jPartIndex*dX[1];
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
              double x = iPartIndex*dX[0];

          
              double position = 0.;
              if(direction == 1) {
                position = x;
              } else if(direction == 2) {
                position = y;
              } else if(direction == 3) {
                position = z;
              }

              double transitionFactor = (1-std::tanh((position-transitionPoint)/transitionWidth))/2.0;
              rhoPtr[bufferIndex]  = rho1+(rho2-rho1)*transitionFactor;
              rhoEPtr[bufferIndex]  = rhoE1+(rhoE2-rhoE1)*transitionFactor;
              for(int iDim=0;iDim<numDim;iDim++){
                size_t velBufferIndex = bufferIndex+numPointsBuffer*iDim;
                if(direction == iDim+1){
                  rhoVPtr[velBufferIndex]  = rho1*velocity1 + 
                    (rho2*velocity2-rho1*velocity1)*transitionFactor;
                }
              }
            }
          }
        }
      } else if (numDim == 2) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];

            double position = 0.;
            if(direction == 1) {
              position = x;
            } else if(direction == 2) {
              position = y;
            }
   
            double transitionFactor = (1-std::tanh((position-transitionPoint)/transitionWidth))/2.0;

            rhoPtr[bufferIndex]  = rho1+(rho2-rho1)*transitionFactor;
            rhoEPtr[bufferIndex]  = rhoE1+(rhoE2-rhoE1)*transitionFactor;
            //if((iIndex == 0 || iIndex == iEnd) && jIndex == jStart) {
              //*messageStream << "rho: " << rhoPtr[bufferIndex] << std::endl;
              //*messageStream << "rhoE internal: " << rhoEPtr[bufferIndex] << std::endl;
            //}
            for(int iDim=0;iDim<numDim;iDim++){
              size_t velBufferIndex = bufferIndex+numPointsBuffer*iDim;
              if(direction == iDim+1){
                rhoVPtr[velBufferIndex]  = rho1*velocity1 + 
                  (rho2*velocity2-rho1*velocity1)*transitionFactor;
                //rhoEPtr[bufferIndex] += rhoVPtr[velBufferIndex]*rhoVPtr[velBufferIndex]/
                  //(2.0*rhoPtr[bufferIndex]);
              }
            }
            //if((iIndex == 0 || iIndex == iEnd) && jIndex == jStart) {
              //*messageStream << "rhoV: " << rhoVPtr[bufferIndex] << std::endl;
              //*messageStream << "rhoE total: " << rhoEPtr[bufferIndex] << std::endl;
            //}

          }
        }
      } else if (numDim == 1) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          rhoPtr[iIndex]  = rho1;
          rhoEPtr[iIndex] = rhoE1;
          rhoVPtr[iIndex] = rhoV[0];
        }
      }
      
      return(0);
    }
    

    template<typename GridType,typename StateType,typename BoundaryType>
    int InitializeProtoY4Test1(const GridType &inGrid,StateType &inState,
                               StateType &inParamState,
                               std::vector<BoundaryType> &inBoundaries,
                               const std::vector<double> &inParams,
                               const std::vector<int> &inFlags,int threadId)
    {

      double pulseAmp       = 0.01;
      double pulseWidth     = 0.000001;
      double pulseRadius    = 0.005;
      double centerX        = 0.5;
      double centerY        = 0.010;
      double crossFlowSpeed = 1.0;
      double injectionSpeed = .2;

      
      int INIT_ACOUSTIC_PULSE    = 1;
      int INIT_CROSSFLOW         = 2;
      int INIT_CROSSFLOW_PROFILE = 4;
      int INIT_INJECTION_FLOW    = 8;
      int INIT_INJECTION_PROFILE = 16;
      int INIT_MASS_GLOB         = 32;
      
      int numParams = inParams.size();
      if(numParams > 0)
        pulseAmp        = inParams[0];
      if(numParams > 1)
        pulseWidth      = inParams[1];
      if(numParams > 2)
        pulseRadius     = inParams[2];
      if(numParams > 3)
        centerX         = inParams[3];
      if(numParams > 4)
        centerY         = inParams[4];
      if(numParams > 5)
        crossFlowSpeed  = inParams[5];
      if(numParams > 6)
        injectionSpeed  = inParams[6];

      int optionBits = 7;
      if(!inFlags.empty()){
        optionBits = inFlags[0];
      }

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      
      

      pcpp::IndexIntervalType bufferInterval;
      bufferInterval.InitSimple(bufferSizes);


      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      std::ostringstream messageStream;
      fixtures::CommunicatorType &gridComm(inGrid.Communicator());

      messageStream << "ProtoY4Test1: Grid Sizes: (";
      pcpp::io::DumpContents(messageStream,gridSizes,",");
      messageStream << ")" << std::endl
                << "ProtoY4Test1: Partition Interval: ";
      partitionInterval.PrettyPrint(messageStream);
      messageStream << std::endl << "ProtoY4Test1: Partition Buffer Interval: ";
      partitionBufferInterval.PrettyPrint(messageStream);
      messageStream << std::endl;

      pcpp::io::Everyone(messageStream.str(),std::cout,gridComm);
      pcpp::io::RenewStream(messageStream);

      std::string pulseType;
      bool pulseOn = false;
      if(optionBits&INIT_ACOUSTIC_PULSE){
        pulseType = "Acoustic";
        pulseOn = true;
      } 
      if(optionBits&INIT_MASS_GLOB){
        pulseType += "Mass";
        pulseOn = true;
      }

      if(!pulseType.empty()){
        messageStream << "ProtoY4Test1: " << pulseType << " pulse amplitude: " 
                  << pulseAmp       << std::endl
                  << "ProtoY4Test1: " 
                  << pulseType << " pulse width:     " 
                  << pulseWidth     << std::endl
                  << "ProtoY4Test1: Relative pulse centerX:         " << centerX        << std::endl
                  << "ProtoY4Test1: Pulse CenterY:         " << centerY        << std::endl;
      } else {
        messageStream << "ProtoY4Test1: No pulse configured." << std::endl;
      }
      if(optionBits&INIT_CROSSFLOW){
        messageStream << "ProtoY4Test1: CrossFlowSpeed:  " << crossFlowSpeed << std::endl;
        if(optionBits&INIT_CROSSFLOW_PROFILE)
          messageStream << "ProtoY4Test1: Crossflow profile enabled." << std::endl;
      }
      if(optionBits&INIT_INJECTION_FLOW){
        messageStream << "ProtoY4Test1: InjectionSpeed:  " << injectionSpeed << std::endl;
        if(optionBits&INIT_INJECTION_PROFILE)
          messageStream << "ProtoY4Test1: Injection flow profile enabled." << std::endl;
      }
      
      if(gridComm.Rank() == 0)
        std::cout << messageStream.str();
      pcpp::io::RenewStream(messageStream);
      
      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      double *gammaPtr      = inParamState.template GetFieldBuffer<double>("gamma");
      double *lengthRefPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *rhoRefPtr     = inParamState.template GetFieldBuffer<double>("refRho");
      double *presRefPtr    = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr    = inParamState.template GetFieldBuffer<double>("refTemperature");
      double *rePtr         = inParamState.template GetFieldBuffer<double>("refRe");
      double *prPtr         = inParamState.template GetFieldBuffer<double>("refPr");
      int *nonDimenPtr      = inParamState.template GetFieldBuffer<int>("nonDimensional");

      double gamma = *gammaPtr;
      double refLength = *lengthRefPtr;
      double refRho = *rhoRefPtr;
      double refPressure = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef = sqrt(gamma*refPressure/refRho);
      double velRef = sqrt(refPressure/refRho);
      if(*nonDimenPtr == 2){
        velRef = sndSpdRef;
        refTemperature *= (gamma-1);
      }
        
      double testSectionP   =   7176.9492187500000;
      double testSectionT   =   121.23677825927734;
      double testSectionRho =   0.22785909473896027;
      double testSectionVel =   566.97581634521487;  // for mach 2.7

      double portP          =   204233.73437500000;
      double portT          =   248.33332824707031;
      double portRho        =   3.1655802726745605;
      double portVel        =   300.53918457031250; // for mach 1.0

      double transitionWidth = 20.0*dX[0];


      // Must have one of these for ProtoY4 configuration
      int combustionChamberIndex = inGrid.GetRegionIndex("ChamberBottom");
      if(combustionChamberIndex < 0)
        combustionChamberIndex = inGrid.GetRegionIndex("Bottom");
      if(combustionChamberIndex < 0)
        return(1);
      
      const simulation::grid::subregion &combustionChamberRegion(gridSubRegions[combustionChamberIndex]);
      const pcpp::IndexIntervalType     &chamberInterval(combustionChamberRegion.regionInterval);
      pcpp::IndexIntervalType combustionChamberInterval(chamberInterval);
      // Augment the boundary extent to create an extent that covers the
      // entire test section.
      combustionChamberInterval[1].first  = gridInterval[1].first;
      combustionChamberInterval[1].second = gridInterval[1].second;

      pcpp::IndexIntervalType chamberPartitionInterval;
      partitionInterval.Overlap(combustionChamberInterval,chamberPartitionInterval);

      if(chamberPartitionInterval.NNodes() > 0){ // Points in the chamber
        
        pcpp::IndexIntervalType chamberBufferInterval;
        chamberPartitionInterval.RelativeTranslation(partitionInterval,
                                                     partitionBufferInterval,
                                                     chamberBufferInterval);
        
        double windowWidth   = 5; // leave this many grid spacings on the edge
        
        double chamberStartX = combustionChamberInterval[0].first  * dX[0];
        double chamberEndX   = combustionChamberInterval[0].second * dX[0];
        double chamberWidth  = chamberEndX - chamberStartX;
        double chamberHalf   = chamberWidth/2.0;
        double chamberCenter = chamberStartX + chamberHalf;
        double pulseWindowX1 = chamberStartX + windowWidth*dX[0];
        double pulseWindowX2 = chamberEndX   - windowWidth*dX[0];
        
        std::cout << "TestSection X := [" << chamberStartX << ","
                  << chamberEndX << "]" << std::endl;
        // these only matter for test pulse
        centerX *= chamberWidth;
        centerX += chamberStartX;
        centerX = .000025;

        double pulseWindowWidth = (centerX - pulseWindowX1);
        double pulseWindowY1    = centerY - pulseWindowWidth/2.0;
        double pulseWindowY2    = centerY + pulseWindowWidth/2.0;

        double profileX1        = chamberStartX + windowWidth*dX[0];
        double profileX2        = chamberEndX   - windowWidth*dX[0];
        double profileWidth     = (profileX2 - profileX1)/2.0;

        messageStream << "ProtoY4Test1: Found combustion chamber interval: ";
        combustionChamberInterval.PrettyPrint(messageStream);
        messageStream << ",(" << chamberStartX << "," << chamberEndX << ")" 
                      << std::endl
                      << "ProtoY4Test1: My combustion chamber interval: ";
        chamberPartitionInterval.PrettyPrint(messageStream);
        messageStream << std::endl;
        if(pulseOn){
          messageStream << "ProtoY4Test1: Pulse location = (" << centerX << ","
                        << centerY << ")" << std::endl;
          messageStream << "ProtoY4Test1: Pulse window width = " << pulseWindowWidth
                        << std::endl;
        }

        pcpp::io::RenewStream(messageStream);

        std::vector<size_t> chamberNodes;
        bufferInterval.GetFlatIndices(chamberBufferInterval,chamberNodes);
          
        std::vector<size_t>::iterator nodeIt = chamberNodes.begin();
        while(nodeIt != chamberNodes.end()){
              
          size_t bufferNodeId = *nodeIt++;
          size_t yBufferIndex = bufferNodeId + numPointsBuffer;
          double x = coordinateData[bufferNodeId];
          double y = coordinateData[bufferNodeId+numPointsBuffer];
          double xDist = (x - centerX);
          double yDist = (y - centerY);
          double radDist = std::sqrt(xDist*xDist + yDist*yDist);
          double pulseWeight = Pulse(pulseAmp,pulseWidth,centerX,centerY,0,x,y,0);
          //double velocityProfile = 
          // Parabola(-1.0/(chamberHalf*chamberHalf),chamberCenter,0,0,x,0,0) + 1.0;
          double velocityProfile = Window(transitionWidth,profileWidth,chamberCenter,0,0,x,0,0);
          double pulseProfile    = Window(transitionWidth,pulseWindowWidth,centerX,centerY,0,x,y,0);


          if(!(optionBits&INIT_CROSSFLOW_PROFILE))
            velocityProfile = 1.0;
          if(radDist > pulseRadius)
            pulseWeight = 0.0;

          rhoPtr[bufferNodeId] = testSectionRho/refRho;

          if(optionBits&INIT_MASS_GLOB)
            rhoPtr[bufferNodeId] *= (1.0 + pulseWeight);

          rhoVPtr[bufferNodeId] = 0.0;

          if(optionBits&INIT_CROSSFLOW){
            rhoVPtr[yBufferIndex]  = rhoPtr[bufferNodeId]*testSectionVel*velocityProfile/velRef;
          }

          rhoEPtr[bufferNodeId] = rhoPtr[bufferNodeId]*testSectionT/refTemperature/gamma + 
            0.5*rhoVPtr[yBufferIndex]*rhoVPtr[yBufferIndex]/rhoPtr[bufferNodeId];

          if(optionBits&INIT_ACOUSTIC_PULSE)
            rhoEPtr[bufferNodeId] *= (1.0 + pulseWeight*pulseProfile);
              
        }
      }

      
      int portIndex = inGrid.GetRegionIndex("Port");
      if(portIndex >= 0){
        
        const simulation::grid::subregion &portRegion(gridSubRegions[portIndex]);
        const pcpp::IndexIntervalType &portInterval(portRegion.regionInterval);

        pcpp::IndexIntervalType myPortInterval;
        partitionInterval.Overlap(portInterval,myPortInterval);

        if(myPortInterval.NNodes() > 0){
          
          double portStartY = portInterval[1].first  * dX[1];
          double portEndY   = portInterval[1].second * dX[1];
          double portWidth  = portEndY - portStartY;
          double portHalf   = portWidth/2.0;
          double portCenter = portStartY + portHalf;
          
          messageStream << "ProtoY4Test1: Initializing Port..." << std::endl
                        << "ProtoY4Test1: Port Interval: ";
          portInterval.PrettyPrint(messageStream);
          messageStream << ", Y = [" << portStartY << "," << portEndY << "]" << std::endl;

          pcpp::IndexIntervalType portBufferInterval;
          myPortInterval.RelativeTranslation(partitionInterval,
                                             partitionBufferInterval,
                                             portBufferInterval);
          std::vector<size_t> portNodes;
          bufferInterval.GetFlatIndices(portBufferInterval,portNodes);
          
          std::vector<size_t>::iterator nodeIt = portNodes.begin();
          while(nodeIt != portNodes.end()){


            size_t bufferNodeId = *nodeIt++;
            double y = coordinateData[bufferNodeId+numPointsBuffer];
            //            double injectionProfile  = 
            //             Parabola(-1.0/(portHalf*portHalf),0,portCenter,0,0,y,0) + 1.0;
            double injectionProfile = std::tanh(10000.0*(y - portStartY)) - 
              std::tanh(10000.0*(y-portEndY)) - 1.0;

            rhoPtr[bufferNodeId] =  portRho/refRho;
            rhoVPtr[bufferNodeId+numPointsBuffer] = 0.0;

            if(!(optionBits&INIT_INJECTION_PROFILE))
              injectionProfile = 1.0;

            rhoVPtr[bufferNodeId]    = rhoPtr[bufferNodeId]*portVel/velRef*injectionProfile;

            rhoEPtr[bufferNodeId] = rhoPtr[bufferNodeId]*portT/refTemperature/gamma + 
              0.5*(rhoVPtr[bufferNodeId]*rhoVPtr[bufferNodeId])/rhoPtr[bufferNodeId];
          }
        }

        // Initialize the plug of gas between the port and the test section
        int plugTopId    = inGrid.GetRegionIndex("PortWallTop");
        if(plugTopId < 0){
          std::cout << "ProtoY4Init:Error: Failed to find required boundary [PortWallTop]. Aborting." << std::endl;
          return(1);
        }

        // First, find out if local process owns part of the plug - this procedure
        // is used a lot. Could be a utility.
        const simulation::grid::subregion &portTopRegion(gridSubRegions[plugTopId]);
        const pcpp::IndexIntervalType &portTopInterval(portTopRegion.regionInterval);
        pcpp::IndexIntervalType plugInterval(portInterval);
        plugInterval[0].second = portTopInterval[0].second;
        plugInterval[0].first++;
        pcpp::IndexIntervalType myPlugInterval;
        partitionInterval.Overlap(plugInterval,myPlugInterval);

        if(myPlugInterval.NNodes() > 0){

          pcpp::IndexIntervalType plugBufferInterval;
          myPlugInterval.RelativeTranslation(partitionInterval,
                                             partitionBufferInterval,
                                             plugBufferInterval);
          std::vector<size_t> plugNodes;
          bufferInterval.GetFlatIndices(plugBufferInterval,plugNodes);
          
          std::vector<size_t>::iterator plugNodeIt = plugNodes.begin();
          while(plugNodeIt != plugNodes.end()){
            size_t bufferIndex = *plugNodeIt++;

            rhoPtr[bufferIndex]                  = testSectionRho/refRho;
            rhoVPtr[bufferIndex]                 = 0.0;
            rhoVPtr[bufferIndex+numPointsBuffer] = 0.0;
            rhoEPtr[bufferIndex]                 = rhoPtr[bufferIndex]*testSectionT/refTemperature/gamma;

          }
        } // initialize plug
      } // domain has a port (if not, hrm?)

      return(0);
  };



    template<typename GridType,typename StateType>
    int InitializeConvectingVortex(const GridType &inGrid,StateType &inState,
                                   StateType &inParam,
                                   const std::vector<double> &inParams,
                                   const std::vector<int> &inFlags,int threadId)
    {
      
      double vortexStrength  = 0.01;
      double vortexWidth     = 0.01;
      double vortexCenterX   = 0.5;
      double vortexCenterY   = 0.5;
      double crossFlowX      = 0.5;
      double crossFlowY      = 0.0;


      int INIT_VORTEX            = 1;
      int INIT_CROSSFLOW         = 2;
      int INIT_CROSSFLOW_PROFILE = 4;
      
      std::ostringstream messageStream;
      int numParams = inParams.size();
      if(numParams > 0)
        vortexStrength   = inParams[0];
      if(numParams > 1)
        vortexWidth      = inParams[1];
      if(numParams > 2)
        vortexCenterX    = inParams[2];
      if(numParams > 3)
        vortexCenterY    = inParams[3];
      if(numParams > 4)
        crossFlowX       = inParams[4];
      if(numParams > 5)
        crossFlowY       = inParams[5];

      int optionBits = 7;
      if(!inFlags.empty()){
        optionBits = inFlags[0];
      }

      messageStream << "VortexInit: vortexStrength = " << vortexStrength << std::endl
                << "VortexInit: vortexWidth    = " << vortexWidth    << std::endl
                << "VortexInit: vortexCenterX  = " << vortexCenterX  << std::endl
                << "VortexInit: vortexCenterY  = " << vortexCenterY  << std::endl
                << "VortexInit: crossFlowX     = " << crossFlowX    << std::endl
                << "VortexInit: crossFlowY     = " << crossFlowY    << std::endl
                << "VortexInit: VORTEX is " << (optionBits&INIT_VORTEX ? "ON" : "OFF") << std::endl
                << "VortexInit: CROSSFLOW is " << (optionBits&INIT_CROSSFLOW ? "ON" : "OFF") << std::endl
                << "VortexInit: INJECTION PROFILE is " << (optionBits&INIT_CROSSFLOW_PROFILE ? "ON" : "OFF")
                << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &gridCoordinates(inGrid.CoordinateData());
      double rho = 1.0;

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);

      double L2M1  = 1.0/(vortexWidth*vortexWidth);
      double gamma     = 1.4;


      for(size_t iPoint  = 0;iPoint < numPointsBuffer;iPoint++){
        size_t yIndex    = numPointsBuffer + iPoint;
        double x         = gridCoordinates[iPoint];
        double y         = gridCoordinates[yIndex];
        double r2L       = L2M1*((x-vortexCenterX)*(x-vortexCenterX) 
                                 + (y-vortexCenterY)*(y-vortexCenterY));
        double eFac      = std::exp(.5*(1.0-r2L));
        double vortexVx  = -1.0*vortexStrength*(y-vortexCenterY)*eFac;
        double vortexVy  =  vortexStrength*(x-vortexCenterX)*eFac;

        if(optionBits&INIT_CROSSFLOW){
          vortexVx += crossFlowX;
          vortexVy += crossFlowY;
        }


        double vortexP = (1.0-0.5*(gamma-1.0)*vortexStrength*vortexStrength*eFac*eFac);
        vortexP = std::pow(vortexP,gamma/(gamma-1.0));
        vortexP *= (1.0/gamma);
        double vortexRho = gamma*vortexP;
        vortexRho = std::pow(vortexRho,(1.0/gamma));


        if(!(optionBits&INIT_VORTEX)){
          vortexRho = 1.0;
          vortexVx  = 0.0;
          vortexVy  = 0.0;
          vortexP   = 1.0;
        }

        rhoPtr[iPoint]  = vortexRho;
        rhoVPtr[iPoint] = vortexRho*vortexVx;
        rhoVPtr[yIndex] = vortexRho*vortexVy;
        rhoEPtr[iPoint] = vortexP/(gamma-1.0) 
          + 0.5*vortexRho*(vortexVx*vortexVx+vortexVy*vortexVy);
      }

      return(0);
    }

    template<typename GridType,typename StateType>
    int InitializeAcousticPulse(const GridType &inGrid,StateType &inState,
                                StateType &inParam,const std::vector<double> &inParams,
                                const std::vector<int> &inFlags,int threadId,
                                std::ostream *messageStream)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      double pulseAmp       = 0.01;
      double pulseWidth     = 0.000001;
      double pulseRadius    = 100.0;
      double relCenterX     = 0.5;
      double relCenterY     = 0.5;
      double relCenterZ     = 0.5;
      double crossFlowSpeed = 0.0;
      double injectionSpeed = 0.0;

      
      int INIT_CROSSFLOW         = 1;
      int INIT_CROSSFLOW_PROFILE = 2;
      
      
      int numParams = inParams.size();
      if(numParams > 0)
        pulseAmp        = inParams[0];
      if(numParams > 1)
        pulseWidth      = inParams[1];
      if(numParams > 2)
        pulseRadius     = inParams[2];
      if(numParams > 3)
        relCenterX      = inParams[3];
      if(numParams > 4)
        relCenterY      = inParams[4];
      if(numParams > 5) 
        relCenterZ      = inParams[5];
      if(numParams > 6)
        crossFlowSpeed  = inParams[6];

      int optionBits = 0;
      if(!inFlags.empty()){
        optionBits = inFlags[0];
      }

      size_t numPointsBuffer = inGrid.BufferSize();             
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &coordinateData(inGrid.CoordinateData());
      const std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
      int numDim = gridSizes.size();
      if(dX.size() != numDim){
        *messageStream << "AcousticPulse:Error: Grid not initialized."
                       << std::endl;
        return(1);
      }

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);
      
      std::vector<double> pulseCenter(3,0);
      pulseCenter[0] = relCenterX;
      pulseCenter[1] = relCenterY;
      pulseCenter[2] = relCenterZ;
      for(int iDim = 0;iDim < numDim;iDim++){
        pulseCenter[iDim] *= (physicalExtent[iDim*2+1] - physicalExtent[iDim*2]);
        pulseCenter[iDim] += physicalExtent[iDim*2];
      }
      if(numDim < 3)
        pulseCenter[2] = 0.0;
      if(numDim < 2)
        pulseCenter[1] = 0.0;


      *messageStream << "AcousticPulse: Grid Sizes: (";
      pcpp::io::DumpContents(*messageStream,gridSizes,",");
      *messageStream << ")" << std::endl
                     << "AcousticPulse: Partition Interval: ";
      partitionInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl << "AcousticPulse: Partition Buffer Interval: ";
      partitionBufferInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl;
      
      *messageStream << "AcousticPulse: Amplitude: " 
                     << pulseAmp       << std::endl
                     << "AcousticPulse: Width: " 
                     << pulseWidth     << std::endl
                     << "AcousticPulse: Relative center: (" << relCenterX 
                     << "," << relCenterY << "," << relCenterZ << ")" << std::endl
                     << "AcousticPulse: Center: (" << pulseCenter[0] << ","
                     << pulseCenter[1] << "," << pulseCenter[2] << ")" << std::endl;
      if(optionBits&INIT_CROSSFLOW){
        *messageStream << "AcousticPulse: CrossFlowSpeed:  " << crossFlowSpeed << std::endl;
        if(optionBits&INIT_CROSSFLOW_PROFILE)
          *messageStream << "AcousticPulse: Crossflow profile enabled." << std::endl;
      }

      double *gammaPtr  = inParam.template GetFieldBuffer<double>("gamma");      
      double gammaValue = *gammaPtr;


      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      

      std::vector<double> velocity(numDim,0.0); // really just velocity at this point
      // Comment next line for quiescent,  uncomment for constant flow
      velocity[0] = crossFlowSpeed;
      
      if(numDim == 3){
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        size_t kStart = partitionBufferInterval[2].first;
        size_t kEnd   = partitionBufferInterval[2].second;
        size_t nPlane = bufferSizes[0]*bufferSizes[1];
        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          size_t kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
            size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
              double x = coordinateData[bufferIndex];
              double y = coordinateData[numPointsBuffer+bufferIndex];
              double z = coordinateData[2*numPointsBuffer+bufferIndex];
              double xRad = x - pulseCenter[0];
              double yRad = y - pulseCenter[1];
              double zRad = z - pulseCenter[2];
              double xyzRad = std::sqrt(xRad*xRad+yRad*yRad+zRad*zRad);

              rhoPtr[bufferIndex]  = 1.0;
              rhoVPtr[bufferIndex] = rhoPtr[bufferIndex]*velocity[0];
              rhoVPtr[bufferIndex+numPointsBuffer] = rhoPtr[bufferIndex]*velocity[1];
              rhoVPtr[bufferIndex+2*numPointsBuffer] = rhoPtr[bufferIndex]*velocity[2];
              rhoEPtr[bufferIndex] = (1.0+Pulse(pulseAmp,pulseWidth,
                                                pulseCenter[0],pulseCenter[1],pulseCenter[2],x,y,z))/(gammaValue-1.0);
              rhoEPtr[bufferIndex] += 
                (rhoVPtr[bufferIndex]*rhoVPtr[bufferIndex] +
                 rhoVPtr[bufferIndex+numPointsBuffer]*rhoVPtr[bufferIndex+numPointsBuffer] +
                 rhoVPtr[bufferIndex+2*numPointsBuffer]*rhoVPtr[bufferIndex+2*numPointsBuffer])/
                (2.0*rhoPtr[bufferIndex]);
            }
          }
        }
      } else if (numDim == 2) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart)+partitionInterval[1].first;
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart)+partitionInterval[0].first;
            double x = coordinateData[bufferIndex];
            double y = coordinateData[numPointsBuffer+bufferIndex];
            double xRad = x - pulseCenter[0];
            double yRad = y - pulseCenter[1];
            double xyzRad = std::sqrt(xRad*xRad+yRad*yRad);
            rhoPtr[bufferIndex]  = 1.0;
            rhoVPtr[bufferIndex] = rhoPtr[bufferIndex]*velocity[0];
            rhoVPtr[bufferIndex+numPointsBuffer] = rhoPtr[bufferIndex]*velocity[1];
            rhoEPtr[bufferIndex] = rhoPtr[bufferIndex]*((1.0/(gammaValue*(gammaValue-1.0))) + 
                                                        .5*(velocity[0]*velocity[0] + velocity[1]*velocity[1]));
            rhoEPtr[bufferIndex] = (1.0 + Pulse(pulseAmp,pulseWidth,
                                                pulseCenter[0],pulseCenter[1],.5,x,y,.5))/(gammaValue-1.0);
            rhoEPtr[bufferIndex] +=
              (rhoVPtr[bufferIndex]*rhoVPtr[bufferIndex] +
               rhoVPtr[bufferIndex+numPointsBuffer]*rhoVPtr[bufferIndex+numPointsBuffer])/
              (2.0*rhoPtr[bufferIndex]);
            
          }
        }
      } else if (numDim == 1) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          size_t iPartIndex = (iIndex - iStart)+partitionInterval[0].first;
          double x = coordinateData[iIndex];
          double xRad = x - pulseCenter[0];
          double xyzRad = std::sqrt(xRad*xRad);
          rhoPtr[iIndex]  = 1.0;
          rhoVPtr[iIndex] = rhoPtr[iIndex]*velocity[0];
          rhoEPtr[iIndex] = (1.0 + Pulse(pulseAmp,pulseWidth,pulseCenter[0],.5,.5,x,.5,.5))/(gammaValue-1.0);
          rhoEPtr[iIndex] += (rhoVPtr[iIndex]*rhoVPtr[iIndex])/(2.0*rhoPtr[iIndex]); 
        }
      }

      
      return(0);
    };

    // initialization for scalar advection test problem
    template<typename GridType,typename StateType>
    int InitializeDensityPulse(const GridType &inGrid,StateType &inState, StateType &inParamState,
                             const std::vector<double> &inParams,int threadId,
                             std::ostream *messageStream = NULL)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      *messageStream << "Starting Density Pulse Flow Initialization" << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      int numDim = gridSizes.size();
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());

      // first 2 parameters are pulse amplitude and width
      // force the user to supply an additional 4 or 6 parameters in 2 or 3D to fully
      // define the problem
      int numParams = inParams.size();
      std::vector<double> position(3,0.0); // initial position
      std::vector<double> velocity(3,0.0); // initial velocity
      double pulseAmp=1.0;
      double pulseWidth=1.0;

      if(numParams < 2){
        *messageStream << "Invalid number of parameters, user must" 
                       << "specify an x location and u velocity " << std::endl;
        return(1);
      } else {
        pulseAmp = inParams[0];
        pulseWidth = inParams[1];
      }

      int locationInd = 2;
      int velocityInd = locationInd+numDim;

      if(numDim == 1) {
        if(numParams == 4){
          position[0] = inParams[locationInd];
          velocity[0] = inParams[velocityInd];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an x location and u velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 2) {
        if(numParams == 6){
          position[0] = inParams[locationInd];
          position[1] = inParams[locationInd+1];
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 3) {
        if(numParams == 8){
          position[0] = inParams[locationInd+0];
          position[1] = inParams[locationInd+1];
          position[2] = inParams[locationInd+2];
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
          velocity[2] = inParams[velocityInd+2];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      }

      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParamState.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParamState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr  = inParamState.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr  = inParamState.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *presRefPtr  = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr  = inParamState.template GetFieldBuffer<double>("refTemperature");

      double gamma = *gammaPtr;
      double refLength = *refLengthPtr;
      double refRho = *rhoRefPtr;
      double refPressure = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef = sqrt(gamma*refPressure/refRho);

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      //*messageStream << "Grid Sizes: (";
      //pcpp::io::DumpContents(*messageStream,gridSizes,",");
      //*messageStream << ")" << std::endl
                //<< "Partition Interval: ";
      //partitionInterval.PrettyPrint(*messageStream);
      //*messageStream << std::endl << "Poiseulle: Partition Buffer Interval: ";
      //partitionBufferInterval.PrettyPrint(*messageStream);
      //*messageStream << std::endl;
      //*messageStream << "Poiseulle: pipe height " << height << std::endl;
      //if(numDim > 2)
        //*messageStream << "Poiseulle: pipe width " << width << std::endl;
      //*messageStream << "Poiseulle: pipe length " << length << std::endl;

      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      // flow conditions
      // working in dimensional units
      double rho = refRho;
      double pressure = refPressure; 
      double temperature = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pressure /= refPressure;
        temperature /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pressure /= gamma*refPressure;
        temperature /= (gamma-1)*refTemperature;
      }

      double rhoE = pressure/(gamma-1)/rho;

      *messageStream << "Density Pulse: Initial state" << std::endl;
      *messageStream << "\t nonDimensional: " << nonDimensionalMode << std::endl;
      //*messageStream << "\t Reynolds Number: " << Re << std::endl;
      *messageStream << "\t rho: " << rho << std::endl;
      *messageStream << "\t pressure: " << pressure << std::endl;
      *messageStream << "\t temperature: " << temperature << std::endl;
      *messageStream << "\t rhoE: " << rhoE << std::endl;
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t position[" << iDim << "]: " << position[iDim] << std::endl;
      }
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t u[" << iDim << "]: " << velocity[iDim] << std::endl;
      }

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;
      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        double z = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          z = kPartIndex*dX[2]; 
        }

        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];

            double densityPulse = Pulse(pulseAmp,pulseWidth,
              position[0],position[1],position[2],x,y,z);
            // cutoff for gaussian
            rhoPtr[bufferIndex] = rho;
            if(densityPulse > 1.e-6)
              rhoPtr[bufferIndex] += densityPulse;

            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rhoPtr[bufferIndex]*velocity[iDim]*velocity[iDim];
            }

          }
        }
      }

      return(0);
    };

    // initialization for scalar advection test problem
    template<typename GridType,typename StateType>
    int InitializeGaussianScalar(const GridType &inGrid,StateType &inState, StateType &inParamState,
                             const std::vector<double> &inParams,int threadId,
                             std::ostream *messageStream = NULL)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      *messageStream << "Starting Gaussian Scalar Flow Initialization" << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      int numDim = gridSizes.size();
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());

      // first 2 parameters are pulse amplitude and width
      // force the user to supply an additional 4 or 6 parameters in 2 or 3D to fully
      // define the problem
      int numParams = inParams.size();
      std::vector<double> position(3,0.0); // initial position
      std::vector<double> velocity(3,0.0); // initial velocity
      double pulseAmp=1.0;
      double pulseWidth=1.0;

      if(numParams < 2){
        *messageStream << "Invalid number of parameters, user must" 
                       << "specify an x location and u velocity " << std::endl;
        return(1);
      } else {
        pulseAmp = inParams[0];
        pulseWidth = inParams[1];
      }

      int locationInd = 2;
      int velocityInd = locationInd+numDim;

      if(numDim == 1) {
        if(numParams == 4){
          position[0] = inParams[locationInd];
          velocity[0] = inParams[velocityInd];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an x location and u velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 2) {
        if(numParams == 6){
          position[0] = inParams[locationInd];
          position[1] = inParams[locationInd+1];
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 3) {
        if(numParams == 8){
          position[0] = inParams[locationInd+0];
          position[1] = inParams[locationInd+1];
          position[2] = inParams[locationInd+2];
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
          velocity[2] = inParams[velocityInd+2];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      }

      // setup the pulse offset for multiple species
      // probably could write a loop to do this cleaner
      std::vector< std::vector<double> > pulseOffset;

      //(0,0,0)
      std::vector<double> tempOffset(numDim,0.);
      pulseOffset.push_back(tempOffset);

      //(1,0,0)
      tempOffset[0] = (readGridExtent[1]-readGridExtent[0])-2*position[0];
      pulseOffset.push_back(tempOffset);

      //(0,1,0)
      tempOffset[0] = 0.;
      tempOffset[1] = (readGridExtent[3]-readGridExtent[2])-2*position[1];
      pulseOffset.push_back(tempOffset);

      //(1,1,0)
      tempOffset[0] = (readGridExtent[1]-readGridExtent[0])-2*position[0];
      pulseOffset.push_back(tempOffset);

      if(numDim > 2){
      //(0,0,readGridExtent[2])
        tempOffset[0] = 0.;
        tempOffset[1] = 0.;
        tempOffset[2] = (readGridExtent[5]-readGridExtent[4])-2*position[2];
        pulseOffset.push_back(tempOffset);
  
        //(readGridExtent[0],0,readGridExtent[2])
        tempOffset[0] = (readGridExtent[1]-readGridExtent[0])-2*position[0];
        pulseOffset.push_back(tempOffset);
  
        //(0,readGridExtent[1],readGridExtent[2])
        tempOffset[0] = 0.;
        tempOffset[1] = (readGridExtent[3]-readGridExtent[2])-2*position[1];
        pulseOffset.push_back(tempOffset);
  
        //(readGridExtent[0],readGridExtent[1],readGridExtent[2])
        tempOffset[0] = (readGridExtent[1]-readGridExtent[0])-2*position[0];
        pulseOffset.push_back(tempOffset);
      }


      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParamState.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParamState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr  = inParamState.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr  = inParamState.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *presRefPtr  = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr  = inParamState.template GetFieldBuffer<double>("refTemperature");

      double gamma = *gammaPtr;
      double refLength = *refLengthPtr;
      double refRho = *rhoRefPtr;
      double refPressure = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef = sqrt(gamma*refPressure/refRho);

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      //*messageStream << "Grid Sizes: (";
      //pcpp::io::DumpContents(*messageStream,gridSizes,",");
      //*messageStream << ")" << std::endl
                //<< "Partition Interval: ";
      //partitionInterval.PrettyPrint(*messageStream);
      //*messageStream << std::endl << "Poiseulle: Partition Buffer Interval: ";
      //partitionBufferInterval.PrettyPrint(*messageStream);
      //*messageStream << std::endl;
      //*messageStream << "Poiseulle: pipe height " << height << std::endl;
      //if(numDim > 2)
        //*messageStream << "Poiseulle: pipe width " << width << std::endl;
      //*messageStream << "Poiseulle: pipe length " << length << std::endl;

      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);

      int scalarHandle = inState.GetDataIndex("scalarVars");
      if(scalarHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &scalarData(inState.Field(scalarHandle));
      int numScalars = inState.Meta()[scalarHandle].ncomp;
      double *scalarPtr    = scalarData.Data<double>();
      if(!scalarPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      for(int iScalar = 0;iScalar < numScalars;iScalar++){
        size_t scalarOffset = iScalar*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          scalarPtr[iPoint+scalarOffset] = 0.0;
      }
      
      // flow conditions
      // working in dimensional units
      double rho = refRho;
      double pressure = refPressure; 
      double temperature = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pressure /= refPressure;
        temperature /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pressure /= gamma*refPressure;
        temperature /= (gamma-1)*refTemperature;
      }

      double rhoE = pressure/(gamma-1)/rho;

      *messageStream << "Scalar Advection: Initial state" << std::endl;
      *messageStream << "\t nonDimensional: " << nonDimensionalMode << std::endl;
      //*messageStream << "\t Reynolds Number: " << Re << std::endl;
      *messageStream << "\t rho: " << rho << std::endl;
      *messageStream << "\t pressure: " << pressure << std::endl;
      *messageStream << "\t temperature: " << temperature << std::endl;
      *messageStream << "\t rhoE: " << rhoE << std::endl;
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t position[" << iDim << "]: " << position[iDim] << std::endl;
      }
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t u[" << iDim << "]: " << velocity[iDim] << std::endl;
      }

      for(int iScalar=0;iScalar<numScalars;iScalar++){
        for(int iDim=0;iDim<numDim;iDim++){
          std::cout << "pulseOffset["<<iScalar<<"]["<<iDim<<"] "<<pulseOffset[iScalar][iDim] << std::endl;
        }
      }

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;
      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        double z = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          z = kPartIndex*dX[2]; 
        }

        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];

            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rho*velocity[iDim]*velocity[iDim];
            }

            // when there is more than one scalar, put multiple globs in varying locations
            // the first scalar is treated as background
            if(numScalars > 1) {
              scalarPtr[bufferIndex] = 1.0;
              for(int iScalar=1;iScalar<numScalars;iScalar++){
  
                double scalarPulse = Pulse(pulseAmp,pulseWidth,
                  position[0]+pulseOffset[iScalar-1][0],
                  position[1]+pulseOffset[iScalar-1][1],
                  position[2]+pulseOffset[iScalar-1][2],
                  x,y,z);
                // cutoff for gaussian
                if(scalarPulse > 1.e-6)
                  scalarPtr[bufferIndex+numPointsBuffer*iScalar] += scalarPulse;
              }
            } else {
              scalarPtr[bufferIndex] = 1.0;
  
              double scalarPulse = Pulse(pulseAmp,pulseWidth,
                position[0], position[1], position[2], x,y,z);
              // cutoff for gaussian
              if(scalarPulse > 1.e-6)
                scalarPtr[bufferIndex] += scalarPulse;
            }
          }
        }
      }

      return(0);
    };

    // initialization for scalar advection-diffusion verification
    template<typename GridType,typename StateType>
    int InitializeAdvectionDiffusion(const GridType &inGrid,StateType &inState, StateType &inParamState,
                                     const std::vector<double> &inParams,int threadId,
                                     std::ostream *messageStream = NULL)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      *messageStream << "Advection-Diffusion Initialization" << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t>     &gridSizes(inGrid.GridSizes());
      const std::vector<size_t>     &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());

      int numDim = gridSizes.size();
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());

      // first 2 parameters are pulse amplitude and width
      // force the user to supply an additional 4 or 6 parameters in 2 or 3D to fully
      // define the problem
      int numParams = inParams.size();
      std::vector<double> position(3,0.0); // initial position
      std::vector<double> velocity(3,0.0); // initial velocity
      double waveAmp=.01;
      double waveK=2.0*3.14159265358979323846;

      if(numParams >= 1)
        waveAmp = inParams[0];
      if(numParams >= 2)
        waveK  = inParams[1];
      
      bool hasPosition = numParams >= (2+numDim);
      bool hasVelocity = numParams >= (2+2*numDim);

      int locationInd = 2;
      int velocityInd = locationInd+numDim;

      if(numDim == 1) {
        if(hasPosition)
          position[0] = inParams[locationInd];
        if(hasVelocity)
          velocity[0] = inParams[velocityInd];
      } else if(numDim == 2) {
        if(hasPosition){
          position[0] = inParams[locationInd];
          position[1] = inParams[locationInd+1];
        } 
        if(hasVelocity){
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
        }
      } else if(numDim == 3) {
        if(hasPosition){
          position[0] = inParams[locationInd+0];
          position[1] = inParams[locationInd+1];
          position[2] = inParams[locationInd+2];
        }
        if(hasVelocity){
          velocity[0] = inParams[velocityInd];
          velocity[1] = inParams[velocityInd+1];
          velocity[2] = inParams[velocityInd+2];
        } 
      }


      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParamState.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParamState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr      = inParamState.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr     = inParamState.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *presRefPtr    = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr    = inParamState.template GetFieldBuffer<double>("refTemperature");

      double gamma          = *gammaPtr;
      double refLength      = *refLengthPtr;
      double refRho         = *rhoRefPtr;
      double refPressure    = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef      = sqrt(gamma*refPressure/refRho);

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);

      int scalarHandle = inState.GetDataIndex("scalarVars");
      if(scalarHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &scalarData(inState.Field(scalarHandle));
      int numScalars = inState.Meta()[scalarHandle].ncomp;
      double *scalarPtr    = scalarData.Data<double>();
      if(!scalarPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      for(int iScalar = 0;iScalar < numScalars;iScalar++){
        size_t scalarOffset = iScalar*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          scalarPtr[iPoint+scalarOffset] = 0.0;
      }
      
      // flow conditions
      // working in dimensional units
      double rho         = refRho;
      double pressure    = refPressure; 
      double temperature = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pressure /= refPressure;
        temperature /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pressure /= gamma*refPressure;
        temperature /= (gamma-1)*refTemperature;
      }

      double rhoE = pressure/(gamma-1)/rho;

      *messageStream << "Scalar Advection Diffusion: Initial state" << std::endl
                     << "\t nonDimensional: " << nonDimensionalMode << std::endl
                     << "\t rho: " << rho << std::endl
                     << "\t pressure: " << pressure << std::endl
                     << "\t temperature: " << temperature << std::endl
                     << "\t rhoE: " << rhoE << std::endl;
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t u[" << iDim << "]: " << velocity[iDim] << std::endl;
      }
      *messageStream << "\t Wave Amp: " << waveAmp << std::endl
                     << "\t Wave Number: " << waveK << std::endl;

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;

      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        double z = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          z = kPartIndex*dX[2]; 
        }

        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];
            double waveCoordinate = x;
            if(numDim > 1)
              if(position[1] > 0) waveCoordinate = y;
            if(numDim > 2)
              if(position[2] > 0) waveCoordinate = z;

            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rho*velocity[iDim]*velocity[iDim];
            }
            for(int iScalar = 0;iScalar < numScalars;iScalar++){
              double phase = iScalar * waveK/4.0;
              size_t scalarOffset = iScalar*numPointsBuffer;
              scalarPtr[scalarOffset+bufferIndex] = waveAmp*std::sin(waveK*waveCoordinate+phase);
            }
          }
        }
      }
      return(0);
    };

    // initialization for scalar advection-diffusion verification
    template<typename GridType,typename StateType>
    int InitializeUniformFlow(const GridType &inGrid,
                              StateType &inState, 
                              StateType &inParamState,
                              const std::vector<double> &inParams,int threadId,
                              std::ostream *messageStream = NULL)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      *messageStream << "Uniform Flow Initialization" << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t>     &gridSizes(inGrid.GridSizes());
      const std::vector<size_t>     &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());

      int numDim = gridSizes.size();

      // first 2 parameters are pulse amplitude and width
      // force the user to supply an additional 4 or 6 parameters in 2 or 3D to fully
      // define the problem
      int numParams = inParams.size();
      std::vector<double> velocity(3,0.0); // initial velocity

      for(int iParam = 0;iParam < numDim;iParam++)
        velocity[iParam] = inParams[iParam];

      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParamState.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParamState.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr      = inParamState.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr     = inParamState.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParamState.template GetFieldBuffer<double>("refLength");
      double *presRefPtr    = inParamState.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr    = inParamState.template GetFieldBuffer<double>("refTemperature");

      double gamma          = *gammaPtr;
      double refLength      = *refLengthPtr;
      double refRho         = *rhoRefPtr;
      double refPressure    = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef      = sqrt(gamma*refPressure/refRho);

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0){
        *messageStream << "Required variable (rho) not found in state." << std::endl;
        return(1);
      }
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0){
        *messageStream << "Required variable (rhoV) not found in state." << std::endl;
        return(1);
      }

      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0){
        *messageStream << "Required variable (rhoE) not found in state." << std::endl;
        return(1);
      }
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      // flow conditions
      // working in dimensional units
      double rho         = refRho;
      double pressure    = refPressure; 
      double temperature = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pressure /= refPressure;
        temperature /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pressure /= gamma*refPressure;
        temperature /= (gamma-1)*refTemperature;
      }

      double rhoE = pressure/(gamma-1)/rho;

      *messageStream << "Uniform Flow: Initial state" << std::endl
                     << "\t nonDimensional: " << nonDimensionalMode << std::endl
                     << "\t rho: " << rho << std::endl
                     << "\t pressure: " << pressure << std::endl
                     << "\t temperature: " << temperature << std::endl
                     << "\t rhoE: " << rhoE << std::endl
                     << "\t velocity: (";
      for(int iDim=0;iDim<numDim;iDim++){
        if(iDim > 0) *messageStream << ",";
        *messageStream << velocity[iDim];
      }
      *messageStream << ")" << std::endl;

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;

      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
        }
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rho*velocity[iDim]*velocity[iDim];
            }
          }
        }
      }
      return(0);
    };

    // initialization for scalar advection-diffusion verification
    template<typename DomainType>
    int InitializeAdvectionDiffusion(DomainType &inDomain,
                                     int iGrid,const std::string &caseName,
                                     int threadId,std::ostream *messageStream = NULL)
    {

      typedef typename DomainType::GridType  GridType;
      typedef typename DomainType::StateType StateType;
      fixtures::ConfigurationType &domainConfig(inDomain.Config());
      
      GridType &inGrid(inDomain.Grid(iGrid));
      StateType &inState(inDomain.State(iGrid));
      StateType &inParams(inDomain.Params(iGrid));
      
      std::string initKey(ConfigKey(ConfigKey(inDomain.Name(),"Init"),caseName));
      std::vector<double> initParams(domainConfig.GetValueVector<double>(ConfigKey(initKey,"Params")));
      std::vector<int> initFlags(domainConfig.GetValueVector<int>(ConfigKey(initKey,"Flags")));

      if(!messageStream)
        messageStream = &std::cout;

      *messageStream << "Advection-Diffusion Initialization" << std::endl;

      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t>     &gridSizes(inGrid.GridSizes());
      const std::vector<size_t>     &bufferSizes(inGrid.BufferSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());

      int numDim = gridSizes.size();
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &readGridExtent(inGrid.PhysicalExtent());

      // first 2 parameters are pulse amplitude and width
      // force the user to supply an additional 4 or 6 parameters in 2 or 3D to fully
      // define the problem
      int numParams = initParams.size();
      std::vector<double> position(3,0.0); // initial position
      std::vector<double> velocity(3,0.0); // initial velocity
      double pulseAmp=1.0;
      double pulseWidth=1.0;

      if(numParams < 2){
        *messageStream << "Invalid number of parameters, user must" 
                       << "specify an x location and u velocity " << std::endl;
        return(1);
      } else {
        pulseAmp = initParams[0];
        pulseWidth = initParams[1];
      }

      int locationInd = 2;
      int velocityInd = locationInd+numDim;

      if(numDim == 1) {
        if(numParams == 4){
          position[0] = initParams[locationInd];
          velocity[0] = initParams[velocityInd];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an x location and u velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 2) {
        if(numParams == 6){
          position[0] = initParams[locationInd];
          position[1] = initParams[locationInd+1];
          velocity[0] = initParams[velocityInd];
          velocity[1] = initParams[velocityInd+1];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      } else if(numDim == 3) {
        if(numParams == 8){
          position[0] = initParams[locationInd+0];
          position[1] = initParams[locationInd+1];
          position[2] = initParams[locationInd+2];
          velocity[0] = initParams[velocityInd];
          velocity[1] = initParams[velocityInd+1];
          velocity[2] = initParams[velocityInd+2];
        } else {
          *messageStream << "Invalid number of parameters, user must" 
                         << "specify an <x,y> location and <u,v> velocity " << std::endl;
          return(1);
        }
      }


      // check on dimensionality
      int nonDimensionalMode=0;
      {
        int paramHandle = inParams.GetDataIndex("nonDimensional");
        if(paramHandle >= 0){
          const pcpp::field::dataset::DataBufferType &paramData(inParams.Field(paramHandle));
          const int *paramPtr = paramData.Data<int>();
          nonDimensionalMode = *paramPtr;
        } else {
          nonDimensionalMode = 0;
        }
      }

      // get reference conditions
      double *gammaPtr      = inParams.template GetFieldBuffer<double>("gamma");
      double *rhoRefPtr     = inParams.template GetFieldBuffer<double>("refRho");
      double *refLengthPtr  = inParams.template GetFieldBuffer<double>("refLength");
      double *presRefPtr    = inParams.template GetFieldBuffer<double>("refPressure");
      double *tempRefPtr    = inParams.template GetFieldBuffer<double>("refTemperature");

      double gamma          = *gammaPtr;
      double refLength      = *refLengthPtr;
      double refRho         = *rhoRefPtr;
      double refPressure    = *presRefPtr;
      double refTemperature = *tempRefPtr;
      double sndSpdRef      = sqrt(gamma*refPressure/refRho);

      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);

      const std::vector<double> &coordinateData(inGrid.CoordinateData());

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoEPtr)
        return(1);

      int scalarHandle = inState.GetDataIndex("scalarVars");
      if(scalarHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &scalarData(inState.Field(scalarHandle));
      int numScalars = inState.Meta()[scalarHandle].ncomp;
      double *scalarPtr    = scalarData.Data<double>();
      if(!scalarPtr)
        return(1);

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      for(int iScalar = 0;iScalar < numScalars;iScalar++){
        size_t scalarOffset = iScalar*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          scalarPtr[iPoint+scalarOffset] = 0.0;
      }
      
      // flow conditions
      // working in dimensional units
      double rho         = refRho;
      double pressure    = refPressure; 
      double temperature = refTemperature;

      if(nonDimensionalMode == 1){
        rho /= refRho;
        pressure /= refPressure;
        temperature /= refTemperature;
      } else if (nonDimensionalMode == 2) {
        rho /= refRho;
        pressure /= gamma*refPressure;
        temperature /= (gamma-1)*refTemperature;
      }

      double rhoE = pressure/(gamma-1)/rho;

      *messageStream << "Scalar Advection Diffusion: Initial state" << std::endl;
      *messageStream << "\t nonDimensional: " << nonDimensionalMode << std::endl;
      //*messageStream << "\t Reynolds Number: " << Re << std::endl;
      *messageStream << "\t rho: " << rho << std::endl;
      *messageStream << "\t pressure: " << pressure << std::endl;
      *messageStream << "\t temperature: " << temperature << std::endl;
      *messageStream << "\t rhoE: " << rhoE << std::endl;
      for(int iDim=0;iDim<numDim;iDim++){
        *messageStream << "\t u[" << iDim << "]: " << velocity[iDim] << std::endl;
      }

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;
      size_t kStart = 0;
      size_t kEnd   = 0;
      size_t nPlane = 0;
      if(numDim > 2 ){
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
        nPlane = bufferSizes[0]*bufferSizes[1];
      }

      for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
        size_t kBufferIndex = kIndex*nPlane;
        size_t kPartIndex = 0;
        double z = 0;
        if (numDim > 2){
          kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          z = kPartIndex*dX[2]; 
        }

        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jkBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
            double x = iPartIndex*dX[0];

            rhoPtr[bufferIndex]  = rho;
            rhoEPtr[bufferIndex] = pressure/(gamma-1);
            for (int iDim=0;iDim<numDim;iDim++){
              rhoVPtr[bufferIndex+numPointsBuffer*iDim] = rhoPtr[bufferIndex]*velocity[iDim];
              rhoEPtr[bufferIndex] += 0.5*rho*velocity[iDim]*velocity[iDim];
            }
            scalarPtr[bufferIndex] = .01*std::sin(2*3.14159265358979323846*x);
            
          }
        }
      }
      return(0);
    };

    template<typename GridType,typename StateType>
    int InitializeDensityPulse(const GridType &inGrid,StateType &inState,
                               StateType &inParam,const std::vector<double> &inParams,
                               const std::vector<int> &inFlags,int threadId,
                               std::ostream *messageStream)
    {
      
      if(!messageStream)
        messageStream = &std::cout;

      double pulseAmp       = 0.01;
      double pulseWidth     = 0.000001;
      double relCenterX     = 0.5;
      double relCenterY     = 0.5;
      double relCenterZ     = 0.5;
      double crossFlowSpeed = 0.0;
      double injectionSpeed = 0.0;

      
      int INIT_CROSSFLOW         = 1;
      int INIT_CROSSFLOW_PROFILE = 2;
      
      
      int numParams = inParams.size();
      if(numParams > 0)
        pulseAmp        = inParams[0];
      if(numParams > 1)
        pulseWidth      = inParams[1];
      if(numParams > 2)
        relCenterX      = inParams[2];
      if(numParams > 3)
        relCenterY      = inParams[3];
      if(numParams > 4) 
        relCenterZ      = inParams[4];
      if(numParams > 5)
        crossFlowSpeed  = inParams[5];

      int optionBits = 0;
      if(!inFlags.empty()){
        optionBits = inFlags[0];
      }

      size_t numPointsBuffer = inGrid.BufferSize();             
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &coordinateData(inGrid.CoordinateData());
      const std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
      int numDim = gridSizes.size();


      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);
      
      std::vector<double> pulseCenter(3,0);
      pulseCenter[0] = relCenterX;
      pulseCenter[1] = relCenterY;
      pulseCenter[2] = relCenterZ;
      for(int iDim = 0;iDim < numDim;iDim++){
        pulseCenter[iDim] *= (physicalExtent[iDim*2+1] - physicalExtent[iDim*2]);
        pulseCenter[iDim] += physicalExtent[iDim*2];
      }
      if(numDim < 3)
        pulseCenter[2] = 0.0;
      if(numDim < 2)
        pulseCenter[1] = 0.0;


      *messageStream << "DensityPulse: Grid Sizes: (";
      pcpp::io::DumpContents(*messageStream,gridSizes,",");
      *messageStream << ")" << std::endl
                     << "DensityPulse: Partition Interval: ";
      partitionInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl << "DensityPulse: Partition Buffer Interval: ";
      partitionBufferInterval.PrettyPrint(*messageStream);
      *messageStream << std::endl;
      
      *messageStream << "DensityPulse: Amplitude: " 
                     << pulseAmp       << std::endl
                     << "DensityPulse: Width: " 
                     << pulseWidth     << std::endl
                     << "DensityPulse: Relative center: (" << relCenterX 
                     << "," << relCenterY << "," << relCenterZ << ")" << std::endl
                     << "DensityPulse: Center: (" << pulseCenter[0] << ","
                     << pulseCenter[1] << "," << pulseCenter[2] << ")" << std::endl;
      if(optionBits > 0){
        *messageStream << "DensityPulse: CrossFlowSpeed:  " << crossFlowSpeed << std::endl;
      }

      double *gammaPtr  = inParam.template GetFieldBuffer<double>("gamma");      
      double gammaValue = *gammaPtr;


      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      

      std::vector<double> velocity(numDim,0.0); // really just velocity at this point
      // Comment next line for quiescent,  uncomment for constant flow
      if(optionBits > 0)
        velocity[optionBits-1] = crossFlowSpeed;
      
      if(numDim == 3){

        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        size_t kStart = partitionBufferInterval[2].first;
        size_t kEnd   = partitionBufferInterval[2].second;
        size_t nPlane = bufferSizes[0]*bufferSizes[1];

        for(size_t kIndex = kStart;kIndex <= kEnd;kIndex++){
          size_t kBufferIndex = kIndex*nPlane;
          size_t kPartIndex = (kIndex - kStart) + partitionInterval[2].first;
          double z = kPartIndex*dX[2]; 
          for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
            size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
            size_t jPartIndex = (jIndex - jStart) + partitionInterval[1].first;
            double y = jPartIndex*dX[1];
            for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
              size_t bufferIndex = jkBufferIndex + iIndex;
              size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
              double x = iPartIndex*dX[0];
              double xRad = x - pulseCenter[0];
              double yRad = y - pulseCenter[1];
              double zRad = z - pulseCenter[2];
              double xyzRad = std::sqrt(xRad*xRad+yRad*yRad+zRad*zRad);

              rhoPtr[bufferIndex]  = 1.0+Pulse(pulseAmp,pulseWidth,pulseCenter[0],pulseCenter[1],pulseCenter[2],x,y,z);
              rhoVPtr[bufferIndex] = rhoPtr[bufferIndex]*velocity[0];
              rhoVPtr[bufferIndex+numPointsBuffer] = rhoPtr[bufferIndex]*velocity[1];
              rhoVPtr[bufferIndex+2*numPointsBuffer] = rhoPtr[bufferIndex]*velocity[2];
              rhoEPtr[bufferIndex] = (1.0/(gammaValue-1.0) + .5*(velocity[0]*velocity[0]   +
                                                                 velocity[1]*velocity[1]   +
                                                                 velocity[2]*velocity[2]));
            }
          }
        }
      } else if (numDim == 2) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        size_t jStart = partitionBufferInterval[1].first;
        size_t jEnd   = partitionBufferInterval[1].second;
        for(size_t jIndex = jStart;jIndex <= jEnd;jIndex++){
          size_t jBufferIndex = jIndex*bufferSizes[0];
          size_t jPartIndex = (jIndex - jStart)+partitionInterval[1].first;
          double y = jPartIndex*dX[1];
          for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
            size_t bufferIndex = jBufferIndex + iIndex;
            size_t iPartIndex = (iIndex - iStart)+partitionInterval[0].first;
            double x = iPartIndex*dX[0];
            double xRad = x - pulseCenter[0];
            double yRad = y - pulseCenter[1];
            double xyzRad = std::sqrt(xRad*xRad+yRad*yRad);
            rhoPtr[bufferIndex]  = 1.0 + Pulse(pulseAmp,pulseWidth,pulseCenter[0],pulseCenter[1],.5,x,y,.5);
            rhoVPtr[bufferIndex] = rhoPtr[bufferIndex]*velocity[0];
            rhoVPtr[bufferIndex+numPointsBuffer] = rhoPtr[bufferIndex]*velocity[1];
            rhoEPtr[bufferIndex] = (1.0/(gammaValue-1.0) + 
                                    .5*(velocity[0]*velocity[0] + velocity[1]*velocity[1]));
            
          }
        }
      } else if (numDim == 1) {
        size_t iStart = partitionBufferInterval[0].first;
        size_t iEnd   = partitionBufferInterval[0].second;
        for(size_t iIndex = iStart;iIndex <= iEnd;iIndex++){
          size_t iPartIndex = (iIndex - iStart)+partitionInterval[0].first;
          double x = iPartIndex*dX[0];
          double xRad = x - pulseCenter[0];
          double xyzRad = std::sqrt(xRad*xRad);
          rhoPtr[iIndex]  = 1.0 + Pulse(pulseAmp,pulseWidth,pulseCenter[0],.5,.5,x,.5,.5);
          rhoVPtr[iIndex] = rhoPtr[iIndex]*velocity[0];
          rhoEPtr[iIndex] = (1.0/(gammaValue - 1.0) + .5*velocity[0]*velocity[0]);
        }
      }

      return(0);
    };

    template<typename GridType,typename StateType>
    int InitializeRiemann1D(const GridType &inGrid,StateType &inState,
                            StateType &paramState,const std::vector<double> &inParams,
                            const std::vector<int> &inFlags,int threadId)
    {

      double rhoL   =  0.0;
      double rhoR   =  0.0;
      double uL     =  0.0;
      double uR     =  0.0;
      double pL     =  0.0;
      double pR     =  0.0;
      double x0     = -1.0;
      int jumpDim   = 0;
      
      
      int numParams = inParams.size();
      if(!(numParams >= 7))
        return(1);
      if(numParams > 0)
        rhoL     = inParams[0];
      if(numParams > 1)
        uL       = inParams[1];
      if(numParams > 2)
        pL       = inParams[2];
      if(numParams > 3)
        rhoR     = inParams[3];
      if(numParams > 4)
        uR       = inParams[4];
      if(numParams > 5) 
        pR       = inParams[5];
      if(numParams > 6)
        x0       = inParams[6];

      if(!inFlags.empty()){
        jumpDim = inFlags[0];
      }

      size_t numPointsBuffer = inGrid.BufferSize();             
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &gridCoords(inGrid.CoordinateData());
      const std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
      int numDim = gridSizes.size();


      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);
      

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      

      double *gammaPtr  = paramState.template GetFieldBuffer<double>("gamma");      
      double gamma      = *gammaPtr;

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      
      // "left" state 
      double rhoVL = rhoL*uL;
      double rhoEL = pL/(gamma-1.0) + 0.5*rhoL*uL*uL;

      // "right" state 
      double rhoVR = rhoR*uR;
      double rhoER = pR/(gamma-1.0) + 0.5*rhoR*uR*uR;

      // locations of jumps, using physical coordinates
      double loc1 = physicalExtent[2*(jumpDim-1)]+x0;
      double loc2 = x0;

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;

      // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
      //  computed correctly
      size_t kStart = 0;
      size_t kEnd   = 0;
      if (numDim == 3) {
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
      }

      for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
        size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
        for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
          size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
          for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
            size_t bufferIndex = jBufferIndex + iIndex;

            //*messageStream << "before cloc, i = " << iIndex << ", j = " << jIndex
            //          << ", k = " << kIndex << std::endl;
            //*messageStream << "  bufferIndex = " << bufferIndex << std::endl;
            //*messageStream << "  jumpDim = " << jumpDim << std::endl;
            //*messageStream << "  gridCoords.size() = " << gridCoords.size() << std::endl;
            double cloc = gridCoords[bufferIndex + (jumpDim-1)*numPointsBuffer];
            //*messageStream << "after cloc, i = " << iIndex << ", j = " << jIndex
            //          << ", k = " << kIndex << std::endl;
            
            if (cloc > loc1 && cloc < loc2) {
              rhoPtr[bufferIndex] = rhoL;
              // V is zero in non-jump directions
              //              for (int vDim=0; vDim<numDim; vDim++) {
              rhoVPtr[bufferIndex + (jumpDim-1)*numPointsBuffer] = rhoVL;
              //              }
              rhoEPtr[bufferIndex] = rhoEL;
            } else {
              rhoPtr[bufferIndex] = rhoR;
              // V is zero in non-jump directions
              //              for (int vDim=0; vDim<numDim; vDim++) {
              rhoVPtr[bufferIndex + (jumpDim-1)*numPointsBuffer] = rhoVR;
              //              }
              rhoEPtr[bufferIndex] = rhoER;
            }
          }
        }
      }
      return(0);
    };

    template<typename GridType,typename StateType>
    int InitializeShocktube(const GridType &inGrid,StateType &inState,
                            StateType &paramState,const std::vector<double> &inParams,
                            const std::vector<int> &inFlags,int threadId)
    {

      double M_init = 0.0;

      int numParams = inParams.size();
      if(!(numParams >= 1))
        return(1);
      if(numParams > 0)
        M_init     = inParams[0];

      //if(!inFlags.empty()){
      //  jumpDim = inFlags[0];
      //}

      size_t numPointsBuffer = inGrid.BufferSize();             
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
      const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<double> &gridCoords(inGrid.CoordinateData());
      const std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
      int numDim = gridSizes.size();


      pcpp::IndexIntervalType gridInterval;
      gridInterval.InitSimple(gridSizes);
      

      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      

      double *gammaPtr  = paramState.template GetFieldBuffer<double>("gamma");      
      double gamma      = *gammaPtr;

      // Zero entire buffer(s)
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 0.0;
        rhoEPtr[iPoint] = 0.0;
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      // reference values
      double rhoRef = 0.42869404474479600;
      double TRef = 298.15;
      double R = 286.9;
      double gm1 = gamma - 1;
      double gp1 = gamma + 1;

      // constants for shock smoothing
      double C = 4.5951198501345889;
      double delta1 = 0.05;

      double rho1 = rhoRef;
      double T1 = TRef;
      double p1 = rho1*R*T1;
      double a1 = std::sqrt(gamma*R*T1);
      double M1 = M_init*M_init;

      double rho2 = rhoRef*(gp1*M1)/(gm1*M1 + 2);
      double p2 = p1*(2*gamma*M1 - gm1)/gp1;
      double T2 = T1*(p2*rho1)/(rho2*p1);
      double M2 = (gm1*M1 + 2)/(2*gamma*M1 - gm1);
      double u2 = M_init*a1 - std::sqrt(M2*gamma*R*T2);

      // since we are using a periodic domain, need to define another location
      //  for a second discontinuity
      double xn = -4.5;

      size_t iStart = partitionBufferInterval[0].first;
      size_t iEnd   = partitionBufferInterval[0].second;
      size_t jStart = partitionBufferInterval[1].first;
      size_t jEnd   = partitionBufferInterval[1].second;

      // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
      //  computed correctly
      size_t kStart = 0;
      size_t kEnd   = 0;
      if (numDim == 3) {
        kStart = partitionBufferInterval[2].first;
        kEnd   = partitionBufferInterval[2].second;
      }

      for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
        size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
        for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
          size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
          for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
            size_t bufferIndex = jBufferIndex + iIndex;

            double x = gridCoords[bufferIndex];
            double eta;
            if (x > -0.1) {
              eta = 0.5*(1 + std::tanh(C*x/delta1));
            } else {
              eta = 0.5*(1 + std::tanh(C*(xn-x)/delta1));
            }

            rhoPtr[bufferIndex] = ((1-eta)*rho2 + eta*rho1)/rhoRef;
            // V is zero in y,z directions
            //              for (int vDim=0; vDim<numDim; vDim++) {
            //rhoVPtr[bufferIndex + (vDim-1)*numPointsBuffer] = rhoV;
            rhoVPtr[bufferIndex] = ((1-eta)*rho2*u2)/(rhoRef*a1);
            //              }
            rhoEPtr[bufferIndex] = ((1-eta)*(p2/gm1 + 0.5*rho2*u2*u2) + eta*(p1/gm1))/(rhoRef*a1*a1);
          }
        }
      }
      return(0);
    };

    template<typename GridType,typename StateType,typename BoundaryType>
    int InitializeHoles(const GridType &inGrid,StateType &inState,
                        StateType &inParam,std::vector<BoundaryType> &inBoundaries)
    {
      size_t numPointsBuffer = inGrid.BufferSize();
      const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
      const std::vector<double> &dX(inGrid.DX());
      const std::vector<size_t> &gridSizes(inGrid.GridSizes());
      const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
      const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());

      int numDim = gridSizes.size();

      double *timePtr   = inState.template GetFieldBuffer<double>("simTime");
      double *deltaTPtr = inParam.template GetFieldBuffer<double>("inputDT");
      double *gammaPtr  = inParam.template GetFieldBuffer<double>("gamma");
      double *CFLPtr    = inParam.template GetFieldBuffer<double>("inputCFL");
      double *cRefPtr   = inParam.template GetFieldBuffer<double>("refC");
      
      double cflValue   = *CFLPtr;
      double cRefValue  = *cRefPtr;
      double gammaValue = *gammaPtr;


      int rhoHandle = inState.GetDataIndex("rho");
      if(rhoHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoData(inState.Field(rhoHandle));
      double *rhoPtr    = rhoData.Data<double>();
      if(!rhoPtr)
        return(1);
      
      int rhoVHandle = inState.GetDataIndex("rhoV");
      if(rhoVHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoVData(inState.Field(rhoVHandle));
      double *rhoVPtr    = rhoVData.Data<double>();
      if(!rhoVPtr)
        return(1);
      
      int rhoEHandle = inState.GetDataIndex("rhoE");
      if(rhoEHandle < 0)
        return(1);
      pcpp::field::dataset::DataBufferType &rhoEData(inState.Field(rhoEHandle));
      double *rhoEPtr    = rhoEData.Data<double>();
      if(!rhoPtr)
        return(1);
      

      // Quiescent and non-singular
      for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++){
        rhoPtr[iPoint] = 1.0;
        rhoEPtr[iPoint] = 1.0/(gammaValue-1.0);
      }
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t dimOffset = iDim*numPointsBuffer;
        for(size_t iPoint = 0;iPoint < numPointsBuffer;iPoint++)
          rhoVPtr[iPoint+dimOffset] = 0.0;
      }
      
      return(0);
    };
  }
}
#endif
//       typename std::vector<BoundaryType>::iterator boundIt = inBoundaries.begin();

//       int numBoundaries = 0;
//       while(boundIt != inBoundaries.end()){

//         BoundaryType &domainBoundary(*boundIt++);
//         std::string boundaryName(domainBoundary.Name());
//         int gridRegionID = domainBoundary.GridRegionID();

//         const simulation::grid::subregion &subRegion(gridSubRegions[gridRegionID]);
//         const pcpp::IndexIntervalType &subRegionInterval(subRegion.regionInterval);
//         const pcpp::IndexIntervalType &subRegionPartitionInterval(subRegion.regionPartitionInterval);
//         const pcpp::IndexIntervalType &subRegionPartitionBufferInterval(subRegion.regionPartitionBufferInterval);

//         if(subRegionPartitionInterval.NNodes() == 0){
//           messageStream << "ProtoY4Test1: Empty boundary: " << boundaryName << std::endl;
//           continue;
//         }

        
//         //        if(boundaryName == "HoleLower" || boundaryName == "Bottom" ){
//         if(boundaryName == "ChamberBottom" || boundaryName == "ChamberTop"){

//           numBoundaries++;

//           // Select only "my" part of the combustion chamber
//           //           if(boundaryName == "HoleLower")
//           //             combustionChamberInterval[0].first = subRegionInterval[0].second + 1;
//           combustionChamberInterval[0] = subRegionInterval[0];
//           combustionChamberInterval.Overlap(partitionInterval,chamberPartitionInterval);