State.C

Go to the documentation of this file.

#include "State.H"
#include "PCPPReport.H"

using namespace pcpp::util;

namespace simulation {

  namespace state {
  
    
    int base::ConfigureState(ConfigurationType &inConfig,std::ostream &messageStream)
    {
      std::string configName(inConfig.GetValue("ConfigName"));
      std::ostringstream Ostr;
      if(!configName.empty())
        Ostr << configName << ":";
      Ostr << "Data";
      std::string solutionKey(Ostr.str());
      std::vector<std::string> fieldNames(inConfig.GetValueVector<std::string>(ConfigKey(solutionKey,"Fields")));
      if(fieldNames.empty()){
        messageStream << "No fields specified for domain"
                      << (configName.empty() ? std::string(" ") : std::string(" ("+configName+") "))
                      << "state" << std::endl;
        return(1);
      }
      std::vector<std::string>::iterator fieldIt = fieldNames.begin();
      while(fieldIt != fieldNames.end()){
        std::string fieldName(*fieldIt++);
        std::string fieldKey(solutionKey+std::string(":")+fieldName);
        std::string fieldUnits = inConfig.GetValue<std::string>(ConfigKey(fieldKey,"Units"));
        int numComponents = inConfig.GetValue<int>(ConfigKey(fieldKey,"Components"));
        int dataSize = inConfig.GetValue<int>(ConfigKey(fieldKey,"DataSize"));
        std::string fieldLocation = inConfig.GetValue<std::string>(ConfigKey(fieldKey,"Location"));
        if(dataSize <= 0 || dataSize > 8){
          messageStream << "ConfigureState:Warning: invalid data size for "
                        << fieldName << ", defaulting to double(8)." << std::endl;
    
          dataSize = 8; // default to double data
        }
        if(numComponents <= 0)
          numComponents = 1; // default to scalar field
        if(fieldLocation.empty())
          fieldLocation = "n";
        if(fieldLocation.size() != 1)
          fieldLocation = fieldLocation.substr(0,1);
        fieldLocation[0] = std::tolower(fieldLocation[0]);
        if(fieldLocation != "n" && fieldLocation != "c" &&
           fieldLocation != "m" && fieldLocation != "s" &&
           fieldLocation != "d"){
          messageStream << "ConfigureState:Warning: invalid solution location for "
                        << fieldName << ", defaulting to nodes." << std::endl;
          fieldLocation = "n";
        }
        AddField(fieldName,fieldLocation[0],numComponents,dataSize,fieldUnits);
      }
      return(0);
    };



    void halo::ConfigureHalo(const pcpp::IndexIntervalType &partitionExtent,
                             const std::vector<pcpp::IndexIntervalType> &haloExtents,
                             size_t numHaloVar)
    {
      numStateComponents = numHaloVar;
      globalPartitionExtent = partitionExtent;
      remoteHaloExtents = haloExtents;
    };

    // Simple halo extent creation - takes a size in each (+/-) direction (i.e. 2*NumDim sizes)
    // and creates the extents that specify the remote halo regions 
    std::vector<pcpp::IndexIntervalType>  halo::CreateRemoteHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                        const pcpp::IndexIntervalType &partitionExtent,
                                                                        std::vector<size_t> &haloSizes)
    {
      std::vector<pcpp::IndexIntervalType> haloExtents;
      int numDim = haloSizes.size()/2;
      haloExtents.resize(numDim*2);
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t numNodesDim = partitionExtent[iDim].second - partitionExtent[iDim].first + 1;
        bool leftHaloRequired = false;
        if(haveNeighbor.empty()){ 
          if(partitionExtent[iDim].first > globalExtent[iDim].first)
            leftHaloRequired = true;
        } else if (haveNeighbor[2*iDim]) {
          leftHaloRequired = true;
        }
        if(leftHaloRequired){
          size_t haloSize = haloSizes[iDim*2];
          pcpp::IndexIntervalType leftHaloExtent(partitionExtent);
//        if((partitionExtent[iDim].first - globalExtent[iDim].first) < haloSize) 
//          haloSize = partitionExtent[iDim].first - globalExtent[iDim].first;
          assert(haloSize < numNodesDim);
          for(int hDim = 0;hDim < numDim;hDim++){
            if(hDim == iDim){
              if(partitionExtent[iDim].first == 0){
                leftHaloExtent[hDim].first  = globalExtent[iDim].second - haloSize + 1;
                leftHaloExtent[hDim].second = globalExtent[iDim].second;
              } else {
                leftHaloExtent[hDim].first  = partitionExtent[iDim].first - haloSize;
                leftHaloExtent[hDim].second = leftHaloExtent[hDim].first + haloSize - 1;
              }
            } else {
              leftHaloExtent[hDim] = partitionExtent[hDim];
            }
          }
          haloExtents[iDim*2].Copy(leftHaloExtent);
        }
        bool rightHaloRequired = false;
        if(haveNeighbor.empty()){
          if(partitionExtent[iDim].second < globalExtent[iDim].second)
            rightHaloRequired = true;
        } else if (haveNeighbor[2*iDim+1]) {
          rightHaloRequired = true;
        }
        if(rightHaloRequired){ // right halo will be required
          size_t haloSize = haloSizes[iDim*2+1];
          pcpp::IndexIntervalType rightHaloExtent(partitionExtent);
//        if((globalExtent[iDim].second - partitionExtent[iDim].second) < haloSize) 
//          haloSize = globalExtent[iDim].second - partitionExtent[iDim].second;
          assert(haloSize < numNodesDim);
          for(int hDim = 0;hDim < numDim;hDim++){
            if(hDim == iDim){
              if(partitionExtent[iDim].second == globalExtent[iDim].second){
                rightHaloExtent[hDim].second  = globalExtent[iDim].first + haloSize - 1;
                rightHaloExtent[hDim].first   = globalExtent[iDim].first;
              } else {
                rightHaloExtent[hDim].first  = partitionExtent[iDim].second + 1;
                rightHaloExtent[hDim].second = rightHaloExtent[hDim].first + haloSize - 1;
              }
            } else {
              rightHaloExtent[hDim] = partitionExtent[hDim];
            }
          }
          haloExtents[iDim*2+1].Copy(rightHaloExtent);
        }
      }
      return(haloExtents);
    };

    std::vector<pcpp::IndexIntervalType>  halo::CreateRemoteHaloExtents(std::vector<size_t> &haloSizes)
    {
      std::vector<pcpp::IndexIntervalType> haloExtents(CreateRemoteHaloExtents(globalGridExtent,globalPartitionExtent,haloSizes));
      return(haloExtents);
    }

    // Simple halo extent creation - takes a size in each (+/-) direction (i.e. 2*NumDim sizes)
    // and creates the extents that specify the local halos (i.e. the local data upon which remote procs depend)
    std::vector<pcpp::IndexIntervalType>  halo::CreateLocalHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                       const pcpp::IndexIntervalType &partitionExtent,
                                                                       std::vector<size_t> &haloSizes)
    {
      std::vector<pcpp::IndexIntervalType> haloExtents;
      int numDim = haloSizes.size()/2;
      haloExtents.resize(numDim*2);
      for(int iDim = 0;iDim < numDim;iDim++){
        size_t numNodesDim = partitionExtent[iDim].second - partitionExtent[iDim].first + 1;
        bool leftHaloRequired = false;
        if(haveNeighbor.empty()){ 
          if(partitionExtent[iDim].first > globalExtent[iDim].first) // left halo will be required
            leftHaloRequired = true;
        } else if (haveNeighbor[iDim*2]) {
          leftHaloRequired = true;
        }
        if(leftHaloRequired){
          size_t haloSize = haloSizes[iDim*2];
          pcpp::IndexIntervalType leftHaloExtent(partitionExtent);
//        if((partitionExtent[iDim].first - globalExtent[iDim].first) < haloSize) 
//          haloSize = partitionExtent[iDim].first - globalExtent[iDim].first;
          assert(haloSize < numNodesDim);
          for(int hDim = 0;hDim < numDim;hDim++){
            if(hDim == iDim){
              leftHaloExtent[hDim].first  = partitionExtent[iDim].first;
              leftHaloExtent[hDim].second = leftHaloExtent[hDim].first + haloSize - 1;
            } else {
              leftHaloExtent[hDim] = partitionExtent[hDim];
            }
          }
          haloExtents[iDim*2].Copy(leftHaloExtent);
        }
        bool rightHaloRequired = false;
        if(haveNeighbor.empty()){
          if(partitionExtent[iDim].second < globalExtent[iDim].second) // right halo will be required
            rightHaloRequired = true;
        } else if(haveNeighbor[2*iDim+1]) {
          rightHaloRequired = true;
        }
        if(rightHaloRequired){
          size_t haloSize = haloSizes[iDim*2+1];
          pcpp::IndexIntervalType rightHaloExtent(partitionExtent);
//        if((globalExtent[iDim].second - partitionExtent[iDim].second) < haloSize) 
//          haloSize = globalExtent[iDim].second - partitionExtent[iDim].second;
          assert(haloSize < numNodesDim);
          for(int hDim = 0;hDim < numDim;hDim++){
            if(hDim == iDim){
              rightHaloExtent[hDim].first  = partitionExtent[iDim].second - haloSize + 1;
              rightHaloExtent[hDim].second = partitionExtent[hDim].second;
            } else {
              rightHaloExtent[hDim] = partitionExtent[hDim];
            }
          }
          haloExtents[iDim*2+1].Copy(rightHaloExtent);
        }
      }
      return(haloExtents);
    };

    std::vector<pcpp::IndexIntervalType>  halo::CreateLocalHaloExtents(std::vector<size_t> &haloSizes)
    {
      std::vector<pcpp::IndexIntervalType> haloExtents(CreateLocalHaloExtents(globalGridExtent,globalPartitionExtent,haloSizes));
      return(haloExtents);
    }

    void halo::ConfigureHalo(const pcpp::IndexIntervalType &globalExtent,
                             const pcpp::IndexIntervalType &partitionExtent,
                             const std::vector<pcpp::IndexIntervalType> &haloExtents)
    {
      //      numStateComponents = numHaloVar;
      globalGridExtent = globalExtent;
      globalPartitionExtent = partitionExtent;
      remoteHaloExtents = haloExtents;
    };

    void halo::ConfigureHalo(const pcpp::IndexIntervalType &globalExtent,
                             const pcpp::IndexIntervalType &partitionExtent)
    {
      //      numStateComponents = numHaloVar;
      globalGridExtent = globalExtent;
      globalPartitionExtent = partitionExtent;
    };

    void halo::CreateHaloBuffers()
    {
      if(remoteHaloExtents.empty())
        return;
      DestroyHaloBuffers();
      ownHaloBuffers = true;
      haloBuffers.resize(numStateFields);
      for(int iField = 0;iField < numStateFields;iField++){
        haloBuffers[iField].resize(remoteHaloExtents.size(),NULL);
        std::vector<pcpp::IndexIntervalType>::iterator haloExtentIt = remoteHaloExtents.begin();
        while(haloExtentIt != remoteHaloExtents.end()){
          int haloIndex = haloExtentIt - remoteHaloExtents.begin();
          pcpp::IndexIntervalType &haloExtent(*haloExtentIt++);
          size_t numNodes = haloExtent.NNodes();
          size_t numVar = numStateComponents * numNodes;
          if(numVar > 0){
            haloBuffers[iField][haloIndex] = new double [numVar];
            if(doFill){
              double *haloBufferPtr = haloBuffers[iField][haloIndex];
              for(size_t iVar = 0;iVar < numVar;iVar++)
                haloBufferPtr[iVar] = -100.0;
            }
          }
        }
      }
    };
    
    void halo::DestroySendBuffers(){
      if(ownSendBuffers){
        std::vector<double *>::iterator sbIt = sendBuffers.begin();
        while(sbIt != sendBuffers.end()){
          double *sendBuffer = *sbIt++;
          if(sendBuffer != NULL)
            delete [] sendBuffer;
        }
      }
      sendBuffers.resize(0);
      sendIndices.resize(0);
    };
    
    void halo::DestroyHaloBuffers(){
      if(ownHaloBuffers){
        int numHalos = remoteHaloExtents.size();
        int numHaloBuffers = haloBuffers.size();
        for(int i = 0;i < numHaloBuffers;i++){
          for(int j = 0;j < numHalos;j++)
            if(haloBuffers[i][j]) delete [] haloBuffers[i][j];
          haloBuffers[i].resize(0);
        }
      }
      haloBuffers.resize(0);
    };
    
    void halo::DestroyRecvBuffers(){
      if(ownRecvBuffers){
        std::vector<double *>::iterator sbIt = recvBuffers.begin();
        while(sbIt != recvBuffers.end()){
          double *recvBuffer = *sbIt++;
          if(recvBuffer) delete [] recvBuffer;
        }
      }
      recvBuffers.resize(0);
      //      recvIndices.resize(0);
      haloRecvIndices.resize(0);
    };
  
    int halo::CreateSimpleSendIndices()
    {
      haveSendData = false;
      size_t totalNumSend = 0;
      if(localHaloExtents.empty())
        return(0);

      int numSendIntervals = localHaloExtents.size();
      sendIndices.resize(numSendIntervals);
      for(int i = 0;i < numSendIntervals;i++){
        size_t numPointsSend = localHaloExtents[i].NNodes();
        if(numPointsSend > 0){

          // Local halo extents are wrt the global grid size. We need the flat array of source
          // buffer indices (i.e. wrt the local buffer size) that will be used to populate 
          // the send buffers. If the local buffer includes "halos" or other buffer zones, 
          // then it is a bit more complicated to obtain the source buffer indices.  
          // The following code determines the source buffer indices.
          if(localPartitionExtent.empty()){ // then halos are not integrated, globalPartitionExtent
            // represents the entire local buffer 
            globalPartitionExtent.GetFlatIndices(localHaloExtents[i],sendIndices[i]);

          } else { // the global partition interval is a sub-interval of the local buffer, deal with it
            //assert(!bufferSizes.empty());
            // Adjust to local numbering: subHaloExtent = localHaloExtents - globalPartitionExtent
            pcpp::IndexIntervalType subHaloExtent(localHaloExtents[i] - globalPartitionExtent);
            // Adjust to the buffer numbering: subHaloExtent = subHaloExtent + localPartitionInterval
            subHaloExtent += localPartitionExtent;
            // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
            pcpp::IndexIntervalType bufferInterval;
            bufferInterval.InitSimple(localBufferSizes);
            bufferInterval.GetFlatIndices(subHaloExtent,sendIndices[i]);
          }

          totalNumSend += numPointsSend;
        }
      }
      haveSendData = (totalNumSend > 0);
      return(totalNumSend);
    }

    int halo::CreateSimpleSendBuffers()
    {
      haveSendData = false;
      size_t totalNumSend = 0;
      if(localHaloExtents.empty())
        return(0);
      DestroySendBuffers();
      ownSendBuffers = true;
      int numSendBuffers = localHaloExtents.size();
      sendBuffers.resize(numSendBuffers,NULL);
      sendIndices.resize(numSendBuffers);
      for(int i = 0;i < numSendBuffers;i++){
        size_t numSend = localHaloExtents[i].NNodes() * numStateComponents;
        if(numSend > 0){
          // Local halo extents are wrt the global mesh. We need the flat array of source
          // buffer indices that will be used to populate the send buffers. If the local
          // buffer includes "halos" or other buffer zones, then it is a bit more 
          // complicated to obtain the source buffer indices.  The following code 
          // determines the source buffer indices.
          if(localPartitionExtent.empty()){ // then halos are not integrated, globalPartitionExtent
            // represents the entire local buffer 
            globalPartitionExtent.GetFlatIndices(localHaloExtents[i],sendIndices[i]);
          } else { // the global partition interval is a sub-interval of the local buffer, deal with it
            //assert(!bufferSizes.empty());
            // Adjust to local numbering: subHaloExtent = localHaloExtents - globalPartitionExtent
            pcpp::IndexIntervalType subHaloExtent(localHaloExtents[i] - globalPartitionExtent);
            // Adjust to the buffer numbering: subHaloExtent = subHaloExtent + localPartitionInterval
            subHaloExtent += localPartitionExtent;
            // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
            pcpp::IndexIntervalType bufferInterval;
            bufferInterval.InitSimple(localBufferSizes);
            bufferInterval.GetFlatIndices(subHaloExtent,sendIndices[i]);
          }

          sendBuffers[i] = new double [numSend];
          // std::cout << "SENDADDRESS(" << i << "): " << sendBuffers[i] 
          //        << std::endl;
          if(doFill){
            double *bufPtr = &sendBuffers[i][0];
            for(size_t iSend = 0;iSend < numSend;iSend++)
              bufPtr[iSend] = -100.0*(i+1);
          }
          totalNumSend += numSend;
        }
      }
      haveSendData = (totalNumSend > 0);
      return(totalNumSend);
    }

    int halo::CreateThreadSendIndices(int threadId)
    {
      if(!haveSendData)
        return(0);
      if(localHaloExtents.empty())
        return(0);
      int numSendBuffers = localHaloExtents.size();
      threadSendIndices[threadId].resize(numSendBuffers);
      threadSendBufferIndices[threadId].resize(numSendBuffers);
      pcpp::IndexIntervalType &threadInterval(threadExtents[threadId]);
      for(int i = 0;i < numSendBuffers;i++){
        double *sendBuffer = sendBuffers[i];
        size_t numSendNodes = localHaloExtents[i].NNodes();
        std::vector<size_t> &threadSendInds(threadSendIndices[threadId][i]);
        std::vector<size_t> &threadSendBufferInds(threadSendBufferIndices[threadId][i]);
        if(numSendNodes > 0){

          // 
          // - debugging prints
          //           std::cout << "Num nodes in Halo = " << numSendNodes << std::endl;
          //           std::cout << "Halo Extent: ";
          //           localHaloExtents[i].PrettyPrint(std::cout);
          //           std::cout << std::endl;
          //           std::cout << "Thread Extent: ";
          //           threadInterval.PrettyPrint(std::cout);
          //           std::cout << std::endl;
          // - debugging
          //

          // collide thread Interval with the halo interval; gets the portion of the thread interval
          // that needs to be sent for this particular border
          pcpp::IndexIntervalType threadSendInterval(threadInterval.Overlap(localHaloExtents[i]));
          size_t numThreadSendNodes = threadSendInterval.NNodes();
          if(numThreadSendNodes > 0){
            //
            // debugging
            //             std::cout << "Thread Send Extent: ";
            //             threadSendInterval.PrettyPrint(std::cout);
            //             std::cout << std::endl;
            //             std::cout << "Thread collides with " << numThreadSendNodes << " nodes to send." << std::endl;
            // debugging
            //
 
            // Get the linear thread-specific source indices; these are the indices of the partition
            // buffer from which to source send data for this thread
            if(localPartitionExtent.empty()){ // then partition extent is not a sub-interval of the buffer
              globalPartitionExtent.GetFlatIndices(threadSendInterval,threadSendInds);
            } else { // partition is a subinterval of the local buffer, deal with it
              //assert(!bufferSizes.empty());
              // Adjust to local numbering: extentWRTLocal = extentWRTGlobal - globalPartitionExtent
              pcpp::IndexIntervalType subThreadSendExtent(threadSendInterval - globalPartitionExtent);
              // Adjust to the buffer numbering: extentWRTLocalBuffer = extentWRTLocal + localPartitionInterval
              subThreadSendExtent += localPartitionExtent;
              // Obtain flat indices: bufferExtent.GetFlatIndices(subHaloExtent,sendIndices[i]);
              pcpp::IndexIntervalType bufferInterval;
              bufferInterval.InitSimple(localBufferSizes);
              bufferInterval.GetFlatIndices(subThreadSendExtent,threadSendInds);
            }

            // Get the linear thread-specific target/pack indices; these are the indices of the send buffer
            // to which to source data is copied for this thread
            localHaloExtents[i].GetFlatIndices(threadSendInterval,threadSendBufferInds);
            if(threadSendInds.size() != numThreadSendNodes){
              exit(1);
            }
            // if the collision interval is not empty, touch the send buffer
            // to situate the page in appropriate NUMA region
            size_t *indexPtr = &threadSendBufferInds[0];
            for(int iComponent = 0;iComponent < numStateComponents;iComponent++){
              for(size_t iSend = 0;iSend < numThreadSendNodes;iSend++){
                // Set to some arbitrary value just to "first touch"
                sendBuffer[indexPtr[iSend]+iComponent*numSendNodes] = -100.0*(threadId+1)*(i+1);
              }
            }
          }
        }
      }
      return(0);
    }

  
    int halo::CreateSendBuffers(const std::vector<pcpp::RemoteCollisionType> &inCollisions) 
    { 
      if(inCollisions.empty())
        return(0);
      DestroySendBuffers();
      ownSendBuffers = true;
      int numCollisions = inCollisions.size();
      sendBuffers.resize(numCollisions);
      sendIndices.resize(numCollisions);
      for(int i = 0;i < numCollisions;i++){
        const pcpp::RemoteCollisionType &remoteCollision(inCollisions[i]);
        const pcpp::IndexIntervalType &collisionExtent(remoteCollision.second);
        globalPartitionExtent.GetFlatIndices(collisionExtent,sendIndices[i]);
        //      int numSend = flatIndices.size() * stateVarIndices.size();
        size_t numSend = sendIndices[i].size() * numStateComponents; // numStateVar includes vector field components
        sendBuffers[i] = new double [numSend];
      }
      return(0);
    };
  
  
    int halo::CreateRecvBuffers(const std::vector<pcpp::RemoteCollisionType> &inCollisions) 
    { 
    
      recvCollisions = inCollisions;
    
      if(inCollisions.empty())
        return(0); 
      DestroyRecvBuffers();
      int numRecvBuffers = 0;
    
      // collide the input collision extents with the global halo extents to figure out
      // the required buffer sizes. this allows multiple remote processors to feed the 
      // local halos 
      std::vector<pcpp::RemoteCollisionType>::const_iterator collIt = inCollisions.begin();
      while(collIt != inCollisions.end()){
      
        int collisionIndex = collIt - inCollisions.begin();
        const pcpp::IndexIntervalType &collisionExtent(collIt++->second);
        size_t numRecv = collisionExtent.NNodes();
        std::vector<pcpp::IndexIntervalType>::iterator haloExtentIt = remoteHaloExtents.begin();
      
        while(haloExtentIt != remoteHaloExtents.end()){
        
          int haloIndex = haloExtentIt - remoteHaloExtents.begin();
        
          pcpp::IndexIntervalType &haloExtent(*haloExtentIt++);
          pcpp::IndexIntervalType haloCollision(haloExtent.Overlap(collisionExtent));
        
          int numCollisions = 0;
          std::vector<size_t> recvIndices;
          if(!haloCollision.empty()){
          
            numCollisions++;
            haloExtent.GetFlatIndices(collisionExtent,recvIndices);
            assert(recvIndices.size() == numRecv);
          
            haloRecvIndices.push_back(std::make_pair(haloIndex,recvIndices));
            size_t numRecvItems = numRecv * numStateComponents;
            double *recvBuffer = new double [numRecvItems];
            recvBuffers.push_back(recvBuffer);
          
          }
          assert(numCollisions == 1);
        }
      }
      return(0);
    };

    int halo::CreateSimpleRecvBuffers()
    { 
      haveRecvData = false;
      if(remoteHaloExtents.empty()){
        std::cout << "halo::CreateSimpleRecvBuffers: WARNING No remote halos." << std::endl;
        return(0); 
      }
      DestroyRecvBuffers();
      size_t totalNumRecv = 0;
      int numRecvBuffers = remoteHaloExtents.size();
      recvBuffers.resize(numRecvBuffers,NULL);
      haloRecvIndices.resize(numRecvBuffers);
      for(int i = 0;i < numRecvBuffers;i++){
        std::vector<size_t> recvIndices;
        size_t numRemoteHaloNodes = remoteHaloExtents[i].NNodes();
        if(numRemoteHaloNodes > 0){
          if(localPartitionExtent.empty()){
            remoteHaloExtents[i].GetFlatIndices(remoteHaloExtents[i],recvIndices);
          } else {
            pcpp::IndexIntervalType localBufferExtent;
            localBufferExtent.InitSimple(localBufferSizes);
            localBufferExtent.GetFlatIndices(remoteHaloBufferExtents[i],recvIndices);
          }
          assert(recvIndices.size() == remoteHaloExtents[i].NNodes());
          haloRecvIndices[i].first  = i;
          haloRecvIndices[i].second = recvIndices;
          size_t numRecvItems = remoteHaloExtents[i].NNodes() * numStateComponents;
          if(numRecvItems > 0) {
            recvBuffers[i] = new double [numRecvItems];
            // std::cout << "RECVADDRESS(" << i << "): " << recvBuffers[i] 
            //        << std::endl;
            if(doFill){
              double *recvBufferPtr = recvBuffers[i];
              for(size_t iItem = 0;iItem < numRecvItems;iItem++)
                recvBufferPtr[iItem] = -1.0*(i+1);
            }
          }
          totalNumRecv += numRecvItems;
        } else {
          haloRecvIndices[i].first = -1;
          recvBuffers[i] = NULL;
        }
      }
      haveRecvData = (totalNumRecv > 0);
      return(totalNumRecv);
    };
  
    int halo::CreateSimpleRecvIndices()
    { 
      haveRecvData = false;
      if(remoteHaloExtents.empty()){
        std::cout << "halo::CreateSimpleRecvBuffers: WARNING No remote halos." << std::endl;
        return(0); 
      }
      //      DestroyRecvBuffers();
      size_t totalNumRecv = 0;
      int numRecvIntervals = remoteHaloExtents.size();
      //      recvBuffers.resize(numRecvBuffers,NULL);
      recvIndices.resize(numRecvIntervals);
      for(int i = 0;i < numRecvIntervals;i++){
        std::vector<size_t> &recvIndex(recvIndices[i]);
        size_t numRemoteHaloNodes = remoteHaloExtents[i].NNodes();
        if(numRemoteHaloNodes > 0){
          if(localPartitionExtent.empty()){
            remoteHaloExtents[i].GetFlatIndices(remoteHaloExtents[i],recvIndex);
          } else {
            pcpp::IndexIntervalType localBufferExtent;
            localBufferExtent.InitSimple(localBufferSizes);
            localBufferExtent.GetFlatIndices(remoteHaloBufferExtents[i],recvIndex);
          }
          assert(recvIndex.size() == remoteHaloExtents[i].NNodes());

          totalNumRecv += numRemoteHaloNodes;
        }
      }
      haveRecvData = (totalNumRecv > 0);
      return(totalNumRecv);
    };

    int halo::CreateMessageBuffers(int numComponents)
    {
      int numCommunication = recvIndices.size();
      if(numCommunication != sendIndices.size())
        return(-1);
      int messageId = communicationBuffers.size();
      halo::messagebuffers newMessageBuffers;
      newMessageBuffers.numComponents = numComponents;
      newMessageBuffers.sendBuffers.resize(numCommunication);
      newMessageBuffers.recvBuffers.resize(numCommunication);
      newMessageBuffers.numPointsSend.resize(numCommunication);
      newMessageBuffers.numPointsRecv.resize(numCommunication);
      for(int iComm = 0;iComm < numCommunication;iComm++){
        int numRecv = recvIndices[iComm].size();
        int numSend = sendIndices[iComm].size();
        newMessageBuffers.numPointsSend[iComm] = numSend;
        newMessageBuffers.numPointsRecv[iComm] = numRecv;
        newMessageBuffers.sendBuffers[iComm] = new double [numComponents * numSend];
        newMessageBuffers.recvBuffers[iComm] = new double [numComponents * numRecv];
      }
      communicationBuffers.push_back(newMessageBuffers);
      return(messageId);
    };

    int halo::PackMessageBuffers(const int messageId,const int componentId,
                              const int numComponents,const double *sourceBuffer)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      assert((componentId+numComponents) <= messageBuffers.numComponents);
      int numMessages = sendIndices.size();
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        std::vector<size_t> &sendIndex(sendIndices[iMessage]);
        size_t numSend = sendIndex.size();
        for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
          for(size_t pointIndex = 0;pointIndex < numSend;pointIndex++){
            size_t sendBufferIndex = pointIndex + (componentId+componentIndex)*numSend;
            size_t sourceBufferIndex = componentIndex*numPointsBuffer;
            messageBuffers.sendBuffers[iMessage][sendBufferIndex] = 
              sourceBuffer[sourceBufferIndex+sendIndex[pointIndex]];
          }
        }
      }
      return(0);
    };

    int halo::PackMessageBuffers(const int messageId,const int componentId,
                              const int numComponents,const double *sourceBuffer,
                              const int threadId)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      assert((componentId+numComponents) <= messageBuffers.numComponents);
      int numMessages = sendIndices.size();
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        std::vector<size_t> &threadSendIndex(threadSendIndices[threadId][iMessage]);
        std::vector<size_t> &threadSendBufferIndex(threadSendBufferIndices[threadId][iMessage]);
        std::vector<size_t> &sendIndex(sendIndices[iMessage]);
        size_t numThreadSend = threadSendIndex.size();
        size_t numSend       = sendIndex.size();
        for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
          for(size_t pointIndex = 0;pointIndex < numThreadSend;pointIndex++){
            size_t sendBufferIndex = threadSendBufferIndex[pointIndex] + (componentId+componentIndex)*numSend;
            size_t sourceBufferIndex = componentIndex*numPointsBuffer + threadSendIndex[pointIndex];
            messageBuffers.sendBuffers[iMessage][sendBufferIndex] = sourceBuffer[sourceBufferIndex];
          }
        }
      }
      return(0);
    };

    int halo::SendMessage(const int messageId,const std::vector<int> &neighborRanks,
                          CommunicatorType &inComm)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      int numMessages = neighborRanks.size();
      int numComponents = messageBuffers.numComponents;
      assert(numMessages == messageBuffers.sendBuffers.size());
      int messageTagBase = messageId*100;
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        int numPointsSend  = messageBuffers.numPointsSend[iMessage];
        double *sendBuffer = messageBuffers.sendBuffers[iMessage];
        int numSendItems = numPointsSend*numComponents;
        int remoteRank = neighborRanks[iMessage];
        int messageTag = messageTagBase + iMessage;
        inComm.ASendBuf(sendBuffer,numSendItems,remoteRank,messageTag);        
      }
      return(0);
    };

    int halo::ReceiveMessage(const int messageId,const std::vector<int> &neighborRanks,
                             CommunicatorType &inComm)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      int numMessages = neighborRanks.size();
      int numComponents = messageBuffers.numComponents;
      assert(numMessages == messageBuffers.sendBuffers.size());
      int messageTagBase = messageId*100;
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        int iDir = iMessage%2;
        int numPointsRecv  = messageBuffers.numPointsSend[iMessage];
        double *recvBuffer = messageBuffers.recvBuffers[iMessage];
        int numRecvItems = numPointsRecv*numComponents;
        int remoteRank = neighborRanks[iMessage];
        int messageTag = messageTagBase + iMessage;
        if(iDir == 0)
          messageTag++;
        else
          messageTag--;
        inComm.ARecvBuf(recvBuffer,numRecvItems,remoteRank,messageTag);        
      }
      inComm.WaitAll();
      return(0);
    };

    int halo::UnPackMessageBuffers(const int messageId,const int componentId,
                                   const int numComponents,double *targetBuffer)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      assert((componentId+numComponents) <= messageBuffers.numComponents);
      int numMessages = recvIndices.size();
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        std::vector<size_t> &recvIndex(recvIndices[iMessage]);
        size_t numRecv = recvIndex.size();
        for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
          for(size_t pointIndex = 0;pointIndex < numRecv;pointIndex++){
            size_t recvBufferIndex = pointIndex + (componentId+componentIndex)*numRecv;
            size_t targetBufferIndex = componentIndex*numPointsBuffer + recvIndex[pointIndex];
            targetBuffer[targetBufferIndex] = messageBuffers.recvBuffers[iMessage][recvBufferIndex];
          }
        }
      }
      return(0);
    };

    int halo::UnPackMessageBuffers(const int messageId,const int componentId,
                                   const int numComponents,double *targetBuffer,
                                   const int threadId)
    {
      assert(messageId < communicationBuffers.size());
      halo::messagebuffers &messageBuffers(communicationBuffers[messageId]);
      assert((componentId+numComponents) <= messageBuffers.numComponents);
      int numMessages = recvIndices.size();
      for(int iMessage = 0;iMessage < numMessages;iMessage++){
        std::vector<size_t> &threadRecvIndex(threadHaloBufferIndices[threadId][iMessage]);
        std::vector<size_t> &threadRecvBufferIndex(threadRecvBufferIndices[threadId][iMessage]);
        size_t recvBufferSize = recvIndices[iMessage].size();
        size_t numThreadRecv = threadRecvIndex.size();
        for(int componentIndex = 0;componentIndex < numComponents;componentIndex++){
          for(size_t pointIndex = 0;pointIndex < numThreadRecv;pointIndex++){
            size_t recvBufferIndex = threadRecvBufferIndex[pointIndex] + 
              (componentId+componentIndex)*recvBufferSize;
            size_t targetBufferIndex = componentIndex*numPointsBuffer + threadRecvIndex[pointIndex];
            targetBuffer[targetBufferIndex] = messageBuffers.recvBuffers[iMessage][recvBufferIndex];
          }
        }
      }
      return(0);
    };


    int halo::CreateThreadRecvIndices(int threadId)
    { 
      if(remoteHaloExtents.empty()){
        std::cout << "halo::CreateThreadRecvIndices: WARNING No remote halos." << std::endl;
        return(0); 
      }


      // thread extent (excluding remote halos) wrt global grid size
      pcpp::IndexIntervalType &threadExtent(threadExtents[threadId]);
      // buffer extent includes remote halos
      pcpp::IndexIntervalType &threadBufferExtent(threadBufferExtents[threadId]);


      int numRecvBuffers = remoteHaloExtents.size();
      // if(threadRecvIndices[threadId].empty())
      //   threadRecvIndices[threadId].resize(numRecvBuffers);
      threadRecvBufferIndices[threadId].resize(numRecvBuffers);
      threadHaloBufferIndices[threadId].resize(numRecvBuffers);
//       std::cout << "Thread Extent: ";
//       threadExtent.PrettyPrint(std::cout);
//       std::cout << std::endl;
      for(int i = 0;i < numRecvBuffers;i++){
        // Check to make sure thread needs to recv (i.e. it owns part of part boundary)
        int recvDirection    = i/2;
        int recvOrientation  = i%2;        
        std::vector<size_t> recvIndices;

        pcpp::IndexIntervalType &remoteHaloExtent(remoteHaloExtents[i]);
        pcpp::IndexIntervalType &remoteHaloBufferExtent(remoteHaloBufferExtents[i]);
        
        std::vector<size_t> &threadRecvBufferIndex(threadRecvBufferIndices[threadId][i]);
        std::vector<size_t> &threadHaloBufferIndex(threadHaloBufferIndices[threadId][i]);

        if(!localPartitionExtent.empty()){
          pcpp::IndexIntervalType localBufferExtent;
          localBufferExtent.InitSimple(localBufferSizes);
          pcpp::IndexIntervalType overlapInterval(threadBufferExtent.Overlap(remoteHaloBufferExtent));
          if(overlapInterval.NNodes() > 0){
            remoteHaloBufferExtent.GetFlatIndices(overlapInterval,threadRecvBufferIndex);
            localBufferExtent.GetFlatIndices(overlapInterval,threadHaloBufferIndex);            
          }
        } else {


          size_t numRemoteHaloNodes = remoteHaloExtent.NNodes();
          //         std::cout << "numRemoteHaloNodes = " << numRemoteHaloNodes << std::endl;
          //         std::cout << "HaloExtent(" << i << "): ";
          //         remoteHaloExtent.PrettyPrint(std::cout);
          //         std::cout << std::endl;
          if(numRemoteHaloNodes > 0){
            
            if(recvOrientation){
              if(threadExtent[recvDirection].second != globalPartitionExtent[recvDirection].second)
                continue;
            } else {
              if(threadExtent[recvDirection].first != globalPartitionExtent[recvDirection].first)
                continue;
            }
            
            pcpp::IndexIntervalType checkExtent(threadExtent);
            checkExtent[recvDirection] = remoteHaloExtent[recvDirection];
            pcpp::IndexIntervalType threadRecvExtent(checkExtent.Overlap(remoteHaloExtent));
            size_t totalNumComponents = 0;
            size_t numThreadRecvNodes = threadRecvExtent.NNodes();
            //           std::cout << "numThreadRecvNodes = " << numThreadRecvNodes << std::endl;
            
            if(numThreadRecvNodes > 0){
              
              remoteHaloExtent.GetFlatIndices(threadRecvExtent,threadRecvBufferIndex);
              size_t *threadRecvBufferIndexPtr = &threadRecvBufferIndex[0];
              //            std::vector<size_t> &threadRecvIndex(threadRecvIndices[threadId][i]);
              //            globalPartitionExtent.GetFlatIndices(threadRecvExtent,threadRecvIndex);
              //            size_t totalNumRecvItems = numRemoteHaloNodes * numStateComponents;
              double *recvBufferPtr = recvBuffers[i];
              for(int iField = 0;iField < numStateFields;iField++){
                int numComponents = numFieldComponents[iField];
                double *haloBufferPtr = NULL;
                haloBufferPtr = haloBuffers[iField][i];
                threadHaloBufferIndex = threadRecvBufferIndex;
                if(doFill){
                  for(int iComponent = 0;iComponent < numComponents;iComponent++){
                    size_t componentRecvIndex = numRemoteHaloNodes*totalNumComponents++;
                    size_t componentHaloIndex = iComponent*numRemoteHaloNodes;
                    for(size_t iItem = 0;iItem < numThreadRecvNodes;iItem++){
                      size_t recvIndex = threadRecvBufferIndexPtr[iItem];
                      recvBufferPtr[recvIndex+componentRecvIndex] = -100.0*(threadId+1)+(i+1);
                      if(haloBufferPtr != NULL) // Halo buffers DNE when halos are integrated
                        haloBufferPtr[recvIndex+componentHaloIndex] = -100.0*(threadId+1)+(i+1);
                    }
                  }
                }
              }
            }
          }
        } 
      }
      return(0);
    };
    
    int halo::UnpackReceiveBuffers()
    {
      if(!haveRecvData)
        return(0);

      std::vector<std::pair<int,std::vector<size_t> > >::iterator haloRecvIndicesIt = 
        haloRecvIndices.begin();
    
      while(haloRecvIndicesIt != haloRecvIndices.end()){
      
        int recvIndex = haloRecvIndicesIt - haloRecvIndices.begin();
        std::pair<int,std::vector<size_t> > &haloRecvInfo(*haloRecvIndicesIt++);
        int haloIndex = haloRecvInfo.first; // encodes direction x{0,1},y{2,3},z{4,5}
        if(haloIndex < 0)
          continue;
        std::vector<size_t> &targetIndices(haloRecvInfo.second);
        size_t numRecv = targetIndices.size();
        if(numRecv <= 0)
          continue;
        double *sourceBuffer = recvBuffers[recvIndex];
        size_t numHalo = remoteHaloExtents[haloIndex].NNodes();
        int iVar = 0;
        int recvBufferIndex = 0;
        for(int iField = 0;iField < numStateFields;iField++){
          double *haloBufferPtr = haloBuffers[iField][haloIndex];
          int numComponents = numFieldComponents[iField];
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            for(size_t iRecv = 0;iRecv < numRecv;iRecv++){
              size_t targetIndex = targetIndices[iRecv] + iComponent*numHalo; 
              haloBufferPtr[targetIndex] = sourceBuffer[recvBufferIndex++];
            }
          }
        }
      }
      return(0);
    };

    int halo::UnpackSimpleRecvBuffers()
    {
      if(!haveRecvData){
        std::cout << "halo::UnpackSimpleRecvBuffers: WARNING NO RECEIVE DATA!!" << std::endl;
        return(0);
      }
      
      std::vector<std::pair<int,std::vector<size_t> > >::iterator haloRecvIndicesIt = 
        haloRecvIndices.begin();
      
      while(haloRecvIndicesIt != haloRecvIndices.end()){
        
        //      int recvIndex = haloRecvIndicesIt - haloRecvIndices.begin();
        std::pair<int,std::vector<size_t> > &haloRecvInfo(*haloRecvIndicesIt++);
        int haloIndex = haloRecvInfo.first; // encodes direction x{0,1},y{2,3},z{4,5}
        if(haloIndex < 0)
          continue;
        std::vector<size_t> &targetIndices(haloRecvInfo.second);
        size_t *targetIndPtr = &targetIndices[0];
        size_t numRecv = targetIndices.size();
        double *sourceBuffer = recvBuffers[haloIndex];
        size_t numHalo = remoteHaloExtents[haloIndex].NNodes();
//         std::cout << "Number of halo(" << haloIndex << ") nodes: " 
//                   << numHalo << std::endl; 
//         std::cout << "NumRecv(" << haloIndex << "): " << numRecv << std::endl;
        if(numRecv <= 0)
          continue;
        // std::cout << "halo::UnpackSimpleRecvBuffers: Unpacking recvIndex(" << recvIndex
        //        << ") Direction/recvbufferid(" << haloIndex << ") ADDRESS:" 
        //        << sourceBuffer << std::endl; 
        int recvBufferIndex = 0;
        for(int iField = 0;iField < numStateFields;iField++){
          double *haloBufferPtr = haloBuffers[iField][haloIndex];
          int numComponents = numFieldComponents[iField];
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t componentIndex = iComponent*numHalo;
            for(size_t iRecv = 0;iRecv < numRecv;iRecv++){
              size_t targetIndex = targetIndPtr[iRecv] + componentIndex;
              haloBufferPtr[targetIndex] = sourceBuffer[recvBufferIndex++];
            }
          }
        }
      }
      return(0);
    };
  
    int halo::UnpackSimpleRecvBuffers(int threadId)
    {
      if(!haveRecvData){
        std::cout << "halo::UnpackSimpleRecvBuffers: WARNING NO RECEIVE DATA!!" << std::endl;
        return(0);
      }

      std::vector<std::pair<int,std::vector<size_t> > >::iterator haloRecvIndicesIt = 
        haloRecvIndices.begin();
      
      while(haloRecvIndicesIt != haloRecvIndices.end()){
        
        //      int recvIndex = haloRecvIndicesIt - haloRecvIndices.begin();
        std::pair<int,std::vector<size_t> > &haloRecvInfo(*haloRecvIndicesIt++);
        int haloIndex = haloRecvInfo.first; // encodes direction x{0,1},y{2,3},z{4,5}
        if(haloIndex < 0){
          continue;
        }
        //        std::cout << "HaloIndex/ThreadId: " << haloIndex << "/" << threadId << std::endl;
        std::vector<size_t> &haloIndices(haloRecvInfo.second);
        if(haloIndices.empty()){
          continue;
        }

        std::vector<size_t> &targetIndices(threadRecvBufferIndices[threadId][haloIndex]);
        size_t numRecv = targetIndices.size();
        size_t numHalo = remoteHaloExtents[haloIndex].NNodes();
//         std::cout << "Number of halo(" << haloIndex << "," << threadId << ") nodes: " 
//                   << numHalo << std::endl; 
//         std::cout << "NumRecv(" << haloIndex << "," << threadId 
//                   << "): " << numRecv << std::endl;
        if(numRecv <= 0 || numHalo == 0){
          continue;
        }
        
        size_t *targetIndPtr = &targetIndices[0];
        double *sourceBuffer = recvBuffers[haloIndex];

        // std::cout << "halo::UnpackSimpleRecvBuffers: Unpacking recvIndex(" << recvIndex
        //        << ") Direction/recvbufferid(" << haloIndex << ") ADDRESS:" 
        //        << sourceBuffer << std::endl; 
      
        int componentId = 0;
        for(int iField = 0;iField < numStateFields;iField++){
          double *haloBufferPtr = haloBuffers[iField][haloIndex];
          int numComponents = numFieldComponents[iField];
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t targetComponentIndex = iComponent*numHalo;
            size_t sourceComponentIndex = numHalo*componentId++;
            for(size_t iRecv = 0;iRecv < numRecv;iRecv++){
              size_t targetIndex = targetIndPtr[iRecv] + targetComponentIndex;
              size_t sourceIndex = targetIndPtr[iRecv] + sourceComponentIndex;
              haloBufferPtr[targetIndex] = sourceBuffer[sourceIndex];
            }
          }
        }
      }
      return(0);
    };
  
    int halo::PackSendBuffers(const state::base &inState)
    {
      if(!haveSendData)
        return(0);

      const pcpp::field::dataset::DataContainerType &solnData(inState.Data());
      const pcpp::field::metadataset &solnMetaData(inState.Meta());
      const std::vector<int> &stateFieldIndices(inState.StateFieldIndices());
    
      std::vector<double *>::iterator sendBufIt = sendBuffers.begin();
      std::vector<std::vector<size_t> >::iterator sendIndIt = sendIndices.begin();
    
      while(sendBufIt != sendBuffers.end()){
      
        double *sendBuffer = *sendBufIt++;
        std::vector<size_t> &sendIndex(*sendIndIt++);
        std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();
        size_t iVar = 0;
        size_t numSend = sendIndex.size();
        size_t *sIndex = NULL;
        if(numSend > 0) sIndex = &sendIndex[0];
        
        while(fieldIndexIt != stateFieldIndices.end()){
        
          int fieldIndex = *fieldIndexIt++;
          const pcpp::field::dataset::DataBufferType &fieldData(solnData[fieldIndex]);
          const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);
          const double *sourceBuffer = fieldData.Data<double>();
          size_t sourceSize = fieldData.NItems();
          int numComponents = fieldMetaData.ncomp;
          size_t dataIndex = iVar * numSend;
        
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t sourceIndex = iComponent*sourceSize;
            for(size_t iData = 0;iData < numSend;iData++)
              sendBuffer[dataIndex++] = sourceBuffer[sourceIndex+sIndex[iData]];
            iVar++;
          }
        
        } 
      
      }
      return(0);
    };
  
    int halo::SimpleSend(std::vector<int> &neighborRanks,CommunicatorType &inComm)
    {

//       if(inComm.NProc() <= 1)
//      return(0);

      if(!haveSendData){
        std::cout << "halo::SimpleSend: WARNING: No Send Data!!" << std::endl;
        return(0);
      }

      std::vector<double *> &sendBuffers(this->SendBuffers());
    
      int numSendBuffers = sendBuffers.size();
    
      for(int i = 0;i < numSendBuffers;i++){
        if(sendBuffers[i] == NULL)
          continue;
        pcpp::IndexIntervalType &localHaloExtent(localHaloExtents[i]);
        if(localHaloExtent.empty())
          return(1);
        int remoteRank = neighborRanks[i];
        size_t numPoints  = localHaloExtent.NNodes();
        size_t numSend    = numPoints * this->numStateComponents;
        // std::cout << "Rank(" << inComm.Rank() << ") send in tag[" 
        //        << i << "] direction to rank (" << remoteRank 
        //        << ") ADDRESS: " << sendBuffers[i] << std::endl;
        // if(numSend < 11){
        //   std::cout << "=============" << std::endl;
        //   for(int isend = 0;isend < numSend;isend++){
        //     std::cout << sendBuffers[i][isend] << " ";
        //   }
        //   std::cout << std::endl
        //          << "==============" << std::endl;
        // }
        inComm.ASendBuf(sendBuffers[i],numSend,remoteRank,i); 
     }
      return(0);
    };

//     int halo::PackSimpleSendBuffers(const state::base &inState)
//     {
//       if(!haveSendData) {
//      std::cout << "halo::PackSimpleSendBuffers: WARNING: No Send data!!" 
//                << std::endl;
//      return(0);
//       }

//       const pcpp::field::dataset::DataContainerType &solnData(inState.Data());
//       const pcpp::field::metadataset &solnMetaData(inState.Meta());
//       //      std::vector<double *>::iterator sendBufIt = sendBuffers.begin();
//       std::vector<std::vector<size_t> >::iterator sendIndIt = sendIndices.begin();
//       int numSendBuffers = sendBuffers.size();
//       //      while(sendBufIt != sendBuffers.end()){      
//       for(int iSend = 0;iSend < numSendBuffers;iSend++){
        
//      int bufferIndex = 0;
//      double *sendBuffer = sendBuffers[iSend];

//      std::vector<size_t> &sendIndex(sendIndices[iSend]);

//      std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();

//      int numSend = sendIndex.size();
      
//      while(fieldIndexIt != stateFieldIndices.end()){
        
//        int fieldIndex = *fieldIndexIt++;
//        const pcpp::field::dataset::DataBufferType &fieldData(solnData[fieldIndex]);
//        const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);
//        const double *sourceBuffer = fieldData.Data<double>();

//        int numComponents = fieldMetaData.ncomp;
//        size_t sourceSize = fieldData.NItems()/numComponents;
        
//        for(int iComponent = 0;iComponent < numComponents;iComponent++){
//          int sourceIndex = iComponent*sourceSize;
//          for(int iData = 0;iData < numSend;iData++)
//            sendBuffer[bufferIndex++] = sourceBuffer[sourceIndex+sendIndex[iData]];
//        }
          
//      } 
        
//       }
//       return(0);
//     };

    int halo::PackSimpleSendBuffers(double inA,const state::base &inX,const state::base &inY)
    {
      if(!haveSendData){
        std::cout << "halo::PackSimpleSendBuffersAXPY: WARNING: No send data!!" << std::endl;
        return(0);
      }

      const pcpp::field::dataset::DataContainerType &yData(inY.Data());
      const pcpp::field::dataset::DataContainerType &xData(inX.Data());
      const pcpp::field::metadataset &solnMetaData(inY.Meta());
      std::vector<double *>::iterator sendBufIt = sendBuffers.begin();
      std::vector<std::vector<size_t> >::iterator sendIndIt = sendIndices.begin();

      while(sendBufIt != sendBuffers.end()){      

        size_t bufferIndex = 0;
        double *sendBuffer = *sendBufIt++;
        std::vector<size_t> &sendIndex(*sendIndIt++);
        size_t *sIndex = &sendIndex[0];
        std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();

        size_t numSend = sendIndex.size();
      
        while(fieldIndexIt != stateFieldIndices.end()){
        
          int fieldIndex = *fieldIndexIt++;

          const pcpp::field::dataset::DataBufferType &yFieldData(yData[fieldIndex]);
          const pcpp::field::dataset::DataBufferType &xFieldData(xData[fieldIndex]);
          const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);

          const double *yBuffer = yFieldData.Data<double>();
          const double *xBuffer = xFieldData.Data<double>();

          int numComponents = fieldMetaData.ncomp;
          size_t sourceSize = yFieldData.NItems()/numComponents;
        
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t sourceIndex = iComponent*sourceSize;
            for(size_t iData = 0;iData < numSend;iData++){
              size_t dataIndex = sourceIndex + sIndex[iData];
              sendBuffer[bufferIndex++] = inA*xBuffer[dataIndex]+yBuffer[dataIndex];
            }
          }
        } 
      }
      return(0);
    };
    
    int halo::PackSimpleSendBuffers(const state::base &inX)
    {
      if(!haveSendData){
        std::cout << "halo::PackSimpleSendBuffers2: WARNING: No send data!!" << std::endl;
        return(0);
      }
      
      const pcpp::field::dataset::DataContainerType &xData(inX.Data());
      const pcpp::field::metadataset &solnMetaData(inX.Meta());
      std::vector<double *>::iterator sendBufIt = sendBuffers.begin();
      std::vector<std::vector<size_t> >::iterator sendIndIt = sendIndices.begin();

      while(sendBufIt != sendBuffers.end()){      

        size_t bufferIndex = 0;
        double *sendBuffer = *sendBufIt++;
        std::vector<size_t> &sendIndex(*sendIndIt++);
        size_t *sIndex = &sendIndex[0];
        std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();

        size_t numSend = sendIndex.size();
      
        while(fieldIndexIt != stateFieldIndices.end()){
        
          int fieldIndex = *fieldIndexIt++;

          const pcpp::field::dataset::DataBufferType &xFieldData(xData[fieldIndex]);
          const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);

          const double *xBuffer = xFieldData.Data<double>();

          int numComponents = fieldMetaData.ncomp;
          size_t sourceSize = xFieldData.NItems()/numComponents;
        
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t sourceIndex = iComponent*sourceSize;
            for(size_t iData = 0;iData < numSend;iData++){
              size_t dataIndex = sourceIndex + sIndex[iData];
              sendBuffer[bufferIndex++] = xBuffer[dataIndex];
            }
          }
        } 
      }
      return(0);
    };

    int halo::PackSimpleSendBuffers(const state::base &inX, int threadId)
    {
      if(localHaloExtents.empty())
        return(0);

      if(!haveSendData){
        std::cout << "halo::PackSimpleSendBuffers2: WARNING: No send data!!" << std::endl;
        return(0);
      }
      
      const pcpp::field::dataset::DataContainerType &xData(inX.Data());
      const pcpp::field::metadataset &solnMetaData(inX.Meta());
      std::vector<double *>::iterator sendBufIt = sendBuffers.begin();

      std::vector<std::vector<size_t> >::iterator sendIndIt = threadSendIndices[threadId].begin();
      std::vector<std::vector<size_t> >::iterator sendBufferIndIt = threadSendBufferIndices[threadId].begin();
      std::vector<std::vector<size_t> >::iterator sendAllIndIt = sendIndices.begin();
      
      while(sendBufIt != sendBuffers.end()){      

        double *sendBuffer = *sendBufIt++;
        std::vector<size_t> &sendIndex(*sendIndIt++);
        std::vector<size_t> &sendBufferIndex(*sendBufferIndIt++);
        size_t *sIndex = &sendIndex[0];
        size_t *sBufferIndex = &sendBufferIndex[0];

        int numPackedComponents = 0;
        std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();
        size_t totalNumSend = sendAllIndIt++->size();
        size_t numSend = sendIndex.size();
      
        while(fieldIndexIt != stateFieldIndices.end()){
        
          int fieldIndex = *fieldIndexIt++;

          const pcpp::field::dataset::DataBufferType &xFieldData(xData[fieldIndex]);
          const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);

          const double *xBuffer = xFieldData.Data<double>();

          int numComponents = fieldMetaData.ncomp;
          size_t sourceSize = xFieldData.NItems()/numComponents;
        
          for(int iComponent = 0;iComponent < numComponents;iComponent++){
            size_t sourceIndex = iComponent*sourceSize;
            size_t targetIndex = (numPackedComponents++)*totalNumSend;
            for(size_t iData = 0;iData < numSend;iData++){
              size_t dataIndex = sourceIndex + sIndex[iData];
              size_t bufferIndex = targetIndex + sBufferIndex[iData];
              sendBuffer[bufferIndex] = xBuffer[dataIndex];
            }
          }
        } 
      }
      return(0);
    };
    
//     int halo::PackSimpleSendBuffers(const state::base &inX,int myThreadId)
//     {
//       if(!haveSendData){
//      std::cout << "halo::PackSimpleSendBuffers2: WARNING: No send data!!" << std::endl;
//      return(0);
//       }
      
//       const pcpp::field::dataset::DataContainerType &xData(inX.Data());
//       const pcpp::field::metadataset &solnMetaData(inX.Meta());
//       std::vector<double *>::iterator sendBufIt = sendBuffers.begin();
//       std::vector<std::vector<size_t> >::iterator sendIndIt = sendIndices.begin();

//       while(sendBufIt != sendBuffers.end()){      

//      int bufferIndex = 0;
//      double *sendBuffer = *sendBufIt++;
//         int sendIndex = *sendIndIt++;
//         std::vector<size_t> &threadPackIndex(sendIndex,myThreadId);
//      std::vector<size_t> &threadSendIndex(sendIndex,myThreadId);
//      std::vector<int>::const_iterator fieldIndexIt = stateFieldIndices.begin();

//      int numSend = threadSendIndex.size();
      
//      while(fieldIndexIt != stateFieldIndices.end()){
        
//        int fieldIndex = *fieldIndexIt++;

//        const pcpp::field::dataset::DataBufferType &xFieldData(xData[fieldIndex]);
//        const pcpp::field::metadata &fieldMetaData(solnMetaData[fieldIndex]);

//        const double *xBuffer = xFieldData.Data<double>();

//        int numComponents = fieldMetaData.ncomp;
//        size_t sourceSize = xFieldData.NItems()/numComponents;
        
//        for(int iComponent = 0;iComponent < numComponents;iComponent++){
//          int sourceIndex = iComponent*sourceSize;
//             int bufferIndex = iComponent*numSend;
//          for(int iData = 0;iData < numSend;iData++){
//            int dataIndex = sourceIndex + threadSendIndex[iData];
//               int packIndex = bufferIndex + threadPackIndex[iData];
//            sendBuffer[packIndex] = xBuffer[dataIndex];
//          }
//        }
//      } 
//       }
//       return(0);
//     };

    int halo::Send(CommunicatorType &inComm)
    {
      if(inComm.NProc() <= 1)
        return(0);
      if(!haveSendData)
        return(0);
      std::vector<double *> &sendBuffers(this->SendBuffers());
      std::vector<pcpp::RemoteCollisionType> &sendCollisions(this->SendCollisions());
    
      int numSendBuffers = sendBuffers.size();
    
      for(int i = 0;i < numSendBuffers;i++){
        pcpp::RemoteCollisionType &sendInfo(sendCollisions[i]);
        int remoteRank = sendInfo.first;
        size_t numPoints  = sendInfo.second.NNodes();
        size_t numSend    = numPoints * this->numStateComponents;
        inComm.ASendBuf(sendBuffers[i],numSend,remoteRank);
      }

      return(0);
    };
  
  

    int halo::PostReceives(CommunicatorType &inComm)
    {
      if(inComm.NProc() <= 1)
        return(0);
      if(this->RecvCollisions().empty())
        return(0);
      if(!haveRecvData)
        return(0);
      int numStateComponents = this->numStateComponents;
      std::vector<pcpp::RemoteCollisionType> &recvCollisions(this->RecvCollisions());
      std::vector<double *> &recvBuffers(this->RecvBuffers());
      int numRecvs = recvBuffers.size();
  
      for(int i = 0;i < numRecvs;i++){
        pcpp::RemoteCollisionType &recvInfo(recvCollisions[i]);
        int remoteRank = recvInfo.first;
        size_t numPoints  = recvInfo.second.NNodes();
        size_t numRecv = numPoints * numStateComponents;
        inComm.ARecvBuf(recvBuffers[i],numRecv,remoteRank);
      }

      return(0);
    };

    int halo::PostSimpleReceives(std::vector<int> &sourceRanks,CommunicatorType &inComm)
    {
      //       if(inComm.NProc() <= 1)
      //        return(0);

      if(!haveRecvData){
        std::cout << "halo::PostSimpleRecveives: WARNING: No data to receive!!" << std::endl;
        return(0);
      }

      int numStateComponents = this->numStateComponents;
      std::vector<double *> &recvBuffers(this->RecvBuffers());
      size_t numRecvs = haloRecvIndices.size();

      for(size_t i = 0;i < numRecvs;i++){

        int recvBufferIndex = haloRecvIndices[i].first;

        if(recvBufferIndex < 0){
          continue;
        }

        int remoteRank = sourceRanks[i];

        if(remoteRank < 0){
          if(recvBuffers[recvBufferIndex] != NULL){
            return(1);
          }
          else{
            continue;
          }
        }

        if(recvBuffers[recvBufferIndex] == NULL){
          std::cout << "No receive buffer exists for this neighbor: " << remoteRank << std::endl;
          return(1);
        }
        // std::cout << "Recv buffer index " << (i==recvBufferIndex ? "matches" : "does not match")
        //        << " direction index." << std::endl;
        std::vector<size_t> &recvIndices(haloRecvIndices[i].second);
        size_t numPoints  = recvIndices.size();
        size_t numRecv = numPoints * numStateComponents;
        int recvDir = i%2;
        int tag = 0;
        if(recvDir == 0)
          tag = i + 1;
        if(recvDir == 1)
          tag = i - 1;
        // std::cout << "Rank(" << inComm.Rank() << ") recv local direction [" 
        //        << i << "," << recvDir  << "] in remote direction tag[" << tag << "] from rank (" 
        //        << remoteRank << ") in RecvBuffer ID (" << recvBufferIndex 
        //        << ") ADDRESS: " << recvBuffers[recvBufferIndex] << std::endl;
        inComm.ARecvBuf(recvBuffers[recvBufferIndex],numRecv,remoteRank,tag);
      }
      return(0);
    };

    void halo::ConfigureData(const state::base &inState)
    {
      std::vector<int> selectFields(inState.StateFieldIndices());
      ConfigureData(inState,selectFields);
    };

    void halo::ConfigureData(const state::base &inState,const std::vector<int> &selectFields)
    {
    
      stateFieldIndices = selectFields;
      numStateFields = stateFieldIndices.size();
      numFieldComponents.resize(numStateFields);
      numStateComponents = 0;
      int fieldCount = 0;
      std::vector<int>::iterator fieldIndexIt = stateFieldIndices.begin();
      while(fieldIndexIt != stateFieldIndices.end()){
        int fieldIndex(*fieldIndexIt++);
        numFieldComponents[fieldCount] = inState.Meta()[fieldIndex].ncomp;
        numStateComponents += numFieldComponents[fieldCount];
        fieldCount++;
      }
    };
      
    int halo::CompleteReceives(CommunicatorType &inComm)
    {
      if(!haveRecvData)
        return(0);
      inComm.WaitAll();
      this->UnpackReceiveBuffers();
      return(0);
    };

    int halo::Wait(CommunicatorType &inComm)
    {
      if(!haveRecvData)
        return(0);
      inComm.WaitAll();
      return(0);
    };

    int halo::CompleteSimpleReceives(CommunicatorType &inComm)
    {
      if(!haveRecvData){
        std::cout << "halo::CompleteSimpleReceives: WARNING: No Data!!" 
                  << std::endl;;
        return(0);
      }
      inComm.WaitAll();
      this->UnpackSimpleRecvBuffers();
      return(0);
    };
  
    void halo::ReportSimpleBufferContents(std::ostream &outStream)
    {
      if(!localHaloExtents.empty()){
        int numLocalHalos = localHaloExtents.size();
        for(int iHalo = 0;iHalo < numLocalHalos;iHalo++){
          pcpp::IndexIntervalType &localHaloExtent(localHaloExtents[iHalo]);
          if(!localHaloExtent.empty()){
            size_t nNodesSend = localHaloExtent.NNodes();
            outStream << "Local halo " << iHalo << " send buffer: ";
            if(!(sendBuffers[iHalo] == NULL)){
              outStream << std::endl;
              int tcomp = 0;
              for(int iField = 0;iField < numStateFields;iField++){
                for(int iComp = 0;iComp < numFieldComponents[iField];iComp++){
                  pcpp::report::StructuredBufferContents(outStream,
                                                         sendBuffers[iHalo+(tcomp*nNodesSend)+iComp*nNodesSend],
                                                         localHaloExtent);
                  outStream << std::endl;
                }
                tcomp += numFieldComponents[iField];
              }
            } else {
              outStream << "[NULL]" << std::endl;
            } 
          }
        }
      }
      if(!remoteHaloExtents.empty()){
        int numRemoteHalos = remoteHaloExtents.size();
        for(int iHalo = 0;iHalo < numRemoteHalos;iHalo++){
          pcpp::IndexIntervalType &remoteHaloExtent(remoteHaloExtents[iHalo]);
          if(!remoteHaloExtent.empty()){
            size_t numPointsHalo = remoteHaloExtent.NNodes();
            outStream << "Remote halo " << iHalo << ": ";
            for(int iField = 0;iField < numStateFields;iField++){
              outStream << "Field " << iField << ": ";
              if(!(haloBuffers[iField][iHalo] == NULL)){
                outStream << std::endl;
                for(int iComp = 0;iComp < numFieldComponents[iField];iComp++){
                  pcpp::report::StructuredBufferContents<double>
                    (outStream,haloBuffers[iField][iHalo],remoteHaloExtent);
                  outStream << std::endl;
                }
              } else {
                outStream << "[NULL]" << std::endl;
              }
            }
          }
        }
      }
      if(!remoteHaloExtents.empty()){
        int numRemoteHalos = remoteHaloExtents.size();
        for(int iHalo = 0;iHalo < numRemoteHalos;iHalo++){
          pcpp::IndexIntervalType &remoteHaloExtent(remoteHaloExtents[iHalo]);
          if(!remoteHaloExtent.empty()){
            size_t nNodesRecv = remoteHaloExtent.NNodes();
            outStream << "Remote halo " << iHalo << " recv buffer: ";
            if(!(recvBuffers[iHalo] == NULL)){
              outStream << std::endl;
              int tcomp = 0;
              for(int iField = 0;iField < numStateFields;iField++){
                for(int iComp = 0;iComp < numFieldComponents[iField];iComp++){
                  pcpp::report::StructuredBufferContents(outStream,
                                                         recvBuffers[iHalo+(tcomp*nNodesRecv)+iComp*nNodesRecv],
                                                         remoteHaloExtent);
                  outStream << std::endl;
                }
                tcomp += numFieldComponents[iField];
              }
            } else {
              outStream << "[NULL]" << std::endl;
            } 
          }
        }
      }
    };
    void probe::SetProbeFields(const std::vector<std::string> &inFields)
    {
      if(!fieldIndices.empty())
        fieldIndices.resize(0);
      probeFields = inFields; 
    };
    void probe::SetDelimiter(const std::string &inPointDelimiter)
    {
      pointDelimiter = inPointDelimiter;
    }
    void probe::SetState(const simulation::state::base &inState)
    {
      statePtr = &inState;
      fieldIndices.resize(0);
    };
    void probe::SetDataInterval(const pcpp::IndexIntervalType &inInterval)
    {
      if(!flatProbeIndices.empty())
        flatProbeIndices.resize(0);
      dataInterval = inInterval;
    };
    void probe::SetProbeInterval(const pcpp::IndexIntervalType &inInterval)
    {
      if(!flatProbeIndices.empty())
        flatProbeIndices.resize(0);
      probeInterval = inInterval;
    };
    int probe::Probe(const simulation::state::base &inState,std::ostream &outStream)
    {
      SetState(inState);
      return(Probe(outStream));
    };
    int probe::Probe(std::ostream &outStream)
    {
      if(statePtr == NULL)
        return(1);
      if(fieldIndices.empty()){
        if(InitFieldIndices())
          return(2);
      }
      if(flatProbeIndices.empty()){
        if(InitFlatProbeIndices())
          return(3);
      }
      std::ostringstream Ostr;
      const pcpp::field::dataset::DataContainerType &probeData(statePtr->Data());
      const pcpp::field::metadataset &probeMeta(statePtr->Meta());
      size_t numNodes  = flatProbeIndices.size();
      int numFields = fieldIndices.size();
      for(size_t iNode = 0;iNode < numNodes;iNode++){
        size_t probeIndex = flatProbeIndices[iNode];
        for(int iField = 0;iField < numFields;iField++){
          int fieldIndex = fieldIndices[iField];
          const pcpp::field::metadata &fieldMeta(probeMeta[fieldIndex]);
          const pcpp::field::databuffer &fieldData(probeData[fieldIndex]);
          int numComponents = fieldMeta.ncomp;
          int dsize  = fieldMeta.dsize;
          int nItems = fieldData.NItems();
          int nPerComponent = nItems/numComponents;
          for(int iComp = 0;iComp < numComponents;iComp++){
            if(dsize == 1){
              const char *dataPtr = fieldData.Data<char>();
              outStream << dataPtr[probeIndex+nPerComponent*iComp];
            } else if (dsize == 4){
              const int  *dataPtr = fieldData.Data<int>();
              outStream << dataPtr[probeIndex+nPerComponent*iComp];
            } else if (dsize == 8){
              const double *dataPtr = fieldData.Data<double>();
              outStream << dataPtr[probeIndex+nPerComponent*iComp];
            }
            outStream << " ";
          }
        }
        outStream << pointDelimiter;
      }
      return(0);
    };
    int probe::InitFieldIndices(){
      fieldIndices.resize(0);
      std::vector<std::string>::iterator fieldNameIt = probeFields.begin();
      while(fieldNameIt != probeFields.end()){
        std::string &fieldName(*fieldNameIt++);
        int fieldIndex = statePtr->GetFieldIndex(fieldName);
        if(fieldIndex < 0)
          return(-1);
        fieldIndices.push_back(fieldIndex);
      }
      return(0);
    };
    int probe::InitFlatProbeIndices(){
      flatProbeIndices.resize(0);
      pcpp::IndexIntervalType overlapInterval;
      dataInterval.Overlap(probeInterval,overlapInterval);
      if(overlapInterval != probeInterval)
        return(1);
      dataInterval.GetFlatIndices(probeInterval,flatProbeIndices);
      if(flatProbeIndices.empty())
        return(1);
      numNodes = flatProbeIndices.size();
      return(0);
    };
  };
};