State.H

Go to the documentation of this file.

#ifndef __STATE_H__
#define __STATE_H__

#include "PCPPTypes.H"
#include "PCPPUtil.H"
#include "PCPPFieldData.H"


namespace simulation {

  namespace state {
    
    typedef pcpp::ConfigType ConfigurationType;
    typedef pcpp::CommunicatorType CommunicatorType;
    
    class base : public pcpp::field::dataset 
    {
      
    public:
      
      //    friend simulation::state::base operator*(double,simulation::state::base &);
      
      base() : pcpp::field::dataset(), numStateVar(0), numStateFields(0) {};
      
      base(const base &inBase) : pcpp::field::dataset(inBase), 
                                 numStateVar(inBase.numStateVar),
                                 numStateFields(inBase.numStateFields),
                                 stateFieldIndices(inBase.stateFieldIndices){};
      //                                  domainVarIndices(inBase.domainVarIndices),
      //                                  paramVarIndices(inBase.paramVarIndices){};
      

      inline base &Zero()
      {
        std::vector<int>::iterator indexIt = stateFieldIndices.begin();
        while(indexIt != stateFieldIndices.end()){
          int fieldIndex(*indexIt++);
          DataBufferType &fieldData(data[fieldIndex]);
          size_t numItems          = fieldData.NItems();
          double *srcData = fieldData.Data<double>();
          for(size_t iItem = 0;iItem < numItems;iItem++)
            srcData[iItem] = 0.0;
        }
        return(*this); 
      }
      
      inline base &operator+=(const base &inBase)
      { 
        assert(stateFieldIndices.size() == inBase.stateFieldIndices.size());
        std::vector<int>::iterator indexIt = stateFieldIndices.begin();
        std::vector<int>::const_iterator inIndexIt = inBase.stateFieldIndices.begin();
        while(indexIt != stateFieldIndices.end()){
          int fieldIndex(*indexIt++);
          int inFieldIndex(*inIndexIt++);
          assert(inBase.meta[inFieldIndex] == meta[fieldIndex]);
          DataBufferType &fieldData(data[fieldIndex]);
          size_t numItems          = fieldData.NItems();

          const double *srcData = inBase.data[inFieldIndex].Data<double>();
          double *targData      = data[fieldIndex].Data<double>();
          for(size_t iItem = 0;iItem < numItems;iItem++)
            targData[iItem] += srcData[iItem];
        }
        return(*this); 
      };
      
      inline base &operator-=(const base &inBase)
      {
        
        assert(stateFieldIndices.size() == inBase.stateFieldIndices.size());
        
        std::vector<int>::const_iterator indexIt  = stateFieldIndices.begin();
        std::vector<int>::const_iterator index2It = inBase.stateFieldIndices.begin();
        
        while(indexIt != stateFieldIndices.end()){

          const int myFieldIndex(*indexIt++);
          const int otherFieldIndex(*index2It++);
          assert((meta[myFieldIndex].ncomp == inBase.meta[otherFieldIndex].ncomp) &&
                 (meta[myFieldIndex].dsize == inBase.meta[otherFieldIndex].dsize));

          DataBufferType &fieldData(data[myFieldIndex]);
          const double *srcData    = inBase.data[otherFieldIndex].Data<double>();
          double *targData         = fieldData.Data<double>();

          size_t numItems          = fieldData.NItems();
          for(size_t iItem = 0;iItem < numItems;iItem++){
            targData[iItem] -= srcData[iItem];
          }

        }
        return(*this); 
      };

      inline base operator+(const base &inBase)
      {
        base returnVal;
        returnVal.Copy(*this);
        returnVal += inBase;
        return(returnVal);
      };

      inline base operator-(const base &inBase)
      {
        base returnVal;
        returnVal.Copy(*this);
        returnVal -= inBase;
        return(returnVal);
      };
      
      inline base &operator*=(double inScalar)
      { 
        std::vector<int>::iterator indexIt = stateFieldIndices.begin();
        while(indexIt != stateFieldIndices.end()){
          int fieldIndex(*indexIt++);
          DataBufferType &fieldData(data[fieldIndex]);
          size_t numItems          = fieldData.NItems();
          double *targData      = fieldData.Data<double>();
          for(size_t iItem = 0;iItem < numItems;iItem++)
            targData[iItem] *= inScalar;
        }
        return(*this); 
      };
      
      inline base operator*(double inScalar)
      {
        base returnValue;
        returnValue.Copy(*this);
        returnValue *= inScalar;
        return(returnValue);
      };
      
      inline base &operator=(const base &inState)
      {
        if(this != &inState){
          Copy(inState);
        }
        return(*this);
      };
      
      virtual void Destroy(const std::string &fieldName)
      {
        pcpp::field::dataset::Destroy();
      };
      virtual void Report(std::ostream &outStream)
      {
        bool verbose = false;
        int numField = NumData();
        //        std::cout << "Number of data items in this state object: " << numField << std::endl;
        for(int iField = 0;iField < numField;iField++){

          int dataIndex = iField;

          const std::string &fieldName(meta[dataIndex].name);
          size_t numPoints     = data[dataIndex].NItems();
          size_t numComponents = meta[dataIndex].ncomp;
          size_t numItems      = numPoints/numComponents;
          size_t dsize         = meta[dataIndex].dsize;
          bool stateField      = (std::find(stateFieldIndices.begin(),stateFieldIndices.end(),iField) 
                                  != stateFieldIndices.end());

          outStream << (stateField ? "" : "State ") << "Field: " 
                    << fieldName                                         << std::endl
                    << "number of components: " << numComponents         << std::endl
                    << "number of items:      " << numItems              << std::endl
                    << "number of bytes/item: " << meta[dataIndex].dsize << std::endl;
          
          
          if(dsize == 1){
            signed char *dataBuffer = GetFieldBuffer<signed char>(fieldName);
            outStream << "buffer address:     " << dataBuffer         << std::endl;
            for(int iComponent = 0;iComponent < numComponents;iComponent++){
              std::ostringstream Ostr;
              Ostr << "component-" << iComponent << " ";
              outStream << (numComponents > 1 ? Ostr.str() : "" ) 
                        << "buffer stats:" << std::endl;
              pcpp::util::ReportBufferStats(outStream,numItems,&dataBuffer[iComponent*numItems]);
            }
          } else if (dsize == 4){
            int *dataBuffer = GetFieldBuffer<int>(fieldName);
            outStream << "buffer address:     " << dataBuffer         << std::endl;
            for(int iComponent = 0;iComponent < numComponents;iComponent++){
              std::ostringstream Ostr;
              Ostr << "component-" << iComponent << " ";
              outStream << (numComponents > 1 ? Ostr.str() : "" ) 
                        << "buffer stats:" << std::endl;
              pcpp::util::ReportBufferStats(outStream,numItems,&dataBuffer[iComponent*numItems]);
            }
          } else if (dsize == 8){
            double *dataBuffer = GetFieldBuffer<double>(fieldName);
            outStream << "buffer address:     " << dataBuffer         << std::endl;
            for(int iComponent = 0;iComponent < numComponents;iComponent++){
              std::ostringstream Ostr;
              Ostr << "component-" << iComponent << " ";
              outStream << (numComponents > 1 ? Ostr.str() : "" ) 
                        << "buffer stats:" << std::endl;
              pcpp::util::ReportBufferStats(outStream,numItems,&dataBuffer[iComponent*numItems]);
            }
          }
        }
      };
      
      // virtual pcpp::field::dataset::DataBufferType &Field(const std::string &fieldName)
      // {
      //   return(pcpp::field::dataset::Field(fieldName)); 
      // };
      
      virtual void CopyStateData(const base &inBase)
      {
        Destroy();
        pcpp::field::dataset::Copy(inBase,inBase.stateFieldIndices);
        numStateVar       = inBase.numStateVar;
        numStateFields    = inBase.numStateFields;
        stateFieldIndices.resize(numStateFields,0);
        for(int iField = 0;iField < numStateFields;iField++){
          int inStateFieldIndex = inBase.stateFieldIndices[iField];
          stateFieldIndices[iField] = iField;
        }
        //         domainVarIndices  = inBase.domainVarIndices;
        //         paramVarIndices   = inBase.paramVarIndices;
      };

      virtual void Copy(const base &inBase)
      {
        Destroy();
        pcpp::field::dataset::Copy(inBase);
        numStateVar       = inBase.numStateVar;
        numStateFields    = inBase.numStateFields;
        stateFieldIndices = inBase.stateFieldIndices;
        //         domainVarIndices  = inBase.domainVarIndices;
        //         paramVarIndices   = inBase.paramVarIndices;
      };
      
      virtual void Copy(const base &inBase,const std::string &copyStates)
      {
        Destroy();
        pcpp::field::dataset::Copy(inBase,copyStates);
      };

      virtual void Destroy()
      {
        numStateVar = 0;
        numStateFields = 0;
        stateFieldIndices.resize(0);
        //         domainVarIndices.resize(0);
        //         paramVarIndices.resize(0);
        pcpp::field::dataset::Destroy();
      };
      
      int Alias(const std::string &sourceField,const std::string &fieldAlias)
      {
        int fieldIndex = this->GetFieldIndex(sourceField);
        if(fieldIndex < 0)
          return(1);
        // Copy source meta into new meta
        // Add field/meta to state
        // Set field buffer to source field buffer
        return(1);
      };

      int ConfigureState(ConfigurationType &inConfig,std::ostream &messageStream);
      
      int NumStateFieldComponents(int stateFieldIndex)
      {
        int fieldIndex = stateFieldIndices[stateFieldIndex];
        return(meta[fieldIndex].ncomp);
      };
      
      
      int NumStateFields() { return(numStateFields); };
      
      int NumStateVar() { return(numStateVar); };
      
      std::vector<std::string> GetStateFieldNames()
      {
        std::vector<std::string> stateFieldNames;
        for(int iField = 0;iField < numStateFields;iField++){
          stateFieldNames.push_back(meta[stateFieldIndices[iField]].name);
        }
        return(stateFieldNames);
      };

      inline double *GetStateFieldData(int stateFieldIndex)
      {
        int fieldIndex = stateFieldIndices[stateFieldIndex];
        return(data[fieldIndex].Data<double>());
      };
      
      inline double *GetStateFieldData(const std::string &fieldName)
      {
        int stateFieldIndex = GetStateFieldIndex(fieldName);
        if(stateFieldIndex < 0)
          return((double *)NULL);
        int fieldIndex = stateFieldIndices[stateFieldIndex];
        return(data[fieldIndex].Data<double>());
      };
      
      inline int GetStateFieldIndex(const std::string &fieldName)
      {
        for(int iField = 0;iField < numStateFields;iField++){
          if(meta[stateFieldIndices[iField]].name == fieldName)
            return(iField);
        }
        return(-1);
      };
      
      const std::vector<int> &StateFieldIndices() const
      { return(stateFieldIndices); };
      
      int SetStateFields(const std::vector<int> &inFieldIndices)
      {
        numStateFields = 0;
        stateFieldIndices.resize(0);
        int nVar = 0;
        int numMeta = NumMeta();
        std::vector<int>::const_iterator fieldIndIt = inFieldIndices.begin();
        while(fieldIndIt != inFieldIndices.end()){
          int fieldIndex(*fieldIndIt++);
          if(fieldIndex < numMeta){
            stateFieldIndices.push_back(fieldIndex);
            nVar += meta[fieldIndex].ncomp;
            numStateFields++;
          }
        }
        numStateVar = nVar;
        return(numStateFields);
      };

      int SetStateFields(const std::vector<std::string> &inFieldNames)
      {
        numStateFields = 0;
        stateFieldIndices.resize(0);
        int nVar = 0;
        int numMeta = NumMeta();
        std::vector<std::string>::const_iterator fieldNameIt = inFieldNames.begin();
        while(fieldNameIt != inFieldNames.end()){
          std::string fieldName(*fieldNameIt++);
          int fieldIndex = GetDataIndex(fieldName);
          if(fieldIndex >= 0){
            stateFieldIndices.push_back(fieldIndex);
            nVar += meta[fieldIndex].ncomp;
            numStateFields++;
          }
        }
        numStateVar = nVar;
        return(numStateFields);
      };

      
      int SetStateFields(const std::string &inFieldNames)
      {
        std::vector<std::string> fieldNameVec;
        std::istringstream Istr(inFieldNames);
        std::string fieldName;
        while(Istr >> fieldName){
          if(!fieldName.empty())
            fieldNameVec.push_back(fieldName);
        }
        return(SetStateFields(fieldNameVec));
      };

      std::string ReportMeta()
      {      
        return(Meta().Report());
      };
      
    protected:
      
      int numStateFields;
      int numStateVar; 
      std::vector<int> stateFieldIndices; 
      //       std::vector<int> domainVarIndices;  /// @brief index(id) of domain-global variables
      //       std::vector<int> paramVarIndices;   /// @brief index(id) of parameter (static) data
      
    }; 
    
    int ConfigureMeta(ConfigurationType &inConfig,pcpp::field::metadataset &inMeta,
                      std::ostream &messageStream);
    
    class probe {
    public:
      probe() : statePtr(NULL),pointDelimiter("\n") {};
      probe(pcpp::IndexIntervalType &inInterval) : statePtr(NULL),pointDelimiter("\n"),
                                                   probeInterval(inInterval) {};
      probe(pcpp::IndexIntervalType &inDataInterval,
            pcpp::IndexIntervalType &inProbeInterval) : statePtr(NULL),pointDelimiter("\n"),
                                                        dataInterval(inDataInterval),
                                                        probeInterval(inProbeInterval) {};
      void SetProbeFields(const std::vector<std::string> &inFields);
      void SetDelimiter(const std::string &inPointDelimiter);
      void SetState(const simulation::state::base &inState);
      void SetDataInterval(const pcpp::IndexIntervalType &inInterval);
      void SetProbeInterval(const pcpp::IndexIntervalType &inInterval);
      virtual int Process(){return(0);};
      int Probe(const simulation::state::base &inState,std::ostream &outStream);
      int Probe(std::ostream &outStream);
    protected:
      const simulation::state::base *statePtr;
      std::vector<std::string>       probeFields;
      std::vector<int>               fieldIndices;
      pcpp::IndexIntervalType        dataInterval;
      pcpp::IndexIntervalType        probeInterval;
      std::vector<size_t>            flatProbeIndices;
      std::string                    pointDelimiter;
    private:
      size_t numNodes;
      int InitFieldIndices();
      int InitFlatProbeIndices();
    };

    class halo {
      
    public:

      struct messagebuffers {
        int numComponents;
        std::vector<size_t> numPointsSend;
        std::vector<size_t> numPointsRecv;
        std::vector<double *> sendBuffers;
        std::vector<double *> recvBuffers;
      };

    public:
      int numStateComponents;
      int numStateFields;
      // These indicate the field indices wrt the state that are included
      // in the halo, and the number of components for the given field, resp.
      std::vector<int> stateFieldIndices;
      std::vector<int> numFieldComponents;
      
      halo() : numStateComponents(0), numStateFields(0), ownSendBuffers(false), 
               ownRecvBuffers(false), ownHaloBuffers(false), doFill(false), numThreads(1)
      {};

      halo(const pcpp::IndexIntervalType &inGlobalExtent,const pcpp::IndexIntervalType &inPartitionExtent) :
        numStateComponents(0), numStateFields(0), ownSendBuffers(false),
        ownRecvBuffers(false), ownHaloBuffers(false), doFill(false), globalGridExtent(inGlobalExtent),
        globalPartitionExtent(inPartitionExtent), numThreads(1)
      {};
      
      void SetGridInterval(const pcpp::IndexIntervalType &inExtent) { globalGridExtent = inExtent; };
      void SetPartitionInterval(const pcpp::IndexIntervalType &inExtent) { globalPartitionExtent = inExtent; };

      bool SetFill(bool yesno){ doFill = yesno; return(doFill); };
      bool Fill() {return(doFill); };

      void SetLocalBufferSizes(const std::vector<size_t> &inBufferSizes) 
      { 
        localBufferSizes = inBufferSizes; 
        numPointsBuffer = 1;
        std::vector<size_t>::iterator bsIt = localBufferSizes.begin();
        while(bsIt != localBufferSizes.end())
          numPointsBuffer *= *bsIt++;
      };

      void SetLocalPartitionExtent(const pcpp::IndexIntervalType &inExtent){ localPartitionExtent = inExtent; };

      std::vector<pcpp::IndexIntervalType> CreateRemoteHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                   const pcpp::IndexIntervalType &partitionExtent,
                                                                   std::vector<size_t> &haloSizes);
      std::vector<pcpp::IndexIntervalType> CreateRemoteHaloExtents(std::vector<size_t> &haloSizes);
      
      std::vector<pcpp::IndexIntervalType> CreateLocalHaloExtents(const pcpp::IndexIntervalType &globalExtent,
                                                                  const pcpp::IndexIntervalType &partitionExtent,
                                                                  std::vector<size_t> &haloSizes);
      std::vector<pcpp::IndexIntervalType> CreateLocalHaloExtents(std::vector<size_t> &haloSizes);
      
      // deprecated
      void ConfigureHalo(const pcpp::IndexIntervalType &partitionExtent,
                         const std::vector<pcpp::IndexIntervalType> &haloExtents,
                         size_t numHaloVar);
      void ConfigureHalo(const pcpp::IndexIntervalType &globalExtent,
                         const pcpp::IndexIntervalType &partitionExtent,
                         const std::vector<pcpp::IndexIntervalType> &haloExtents);
      void ConfigureHalo(const pcpp::IndexIntervalType &globalExtent,
                         const pcpp::IndexIntervalType &partitionExtent);
      void ConfigureData(const state::base &inState);
      void ConfigureData(const state::base &inState,const std::vector<int> &selectFields);
      
      // -------- Threaded interface --------------
      int CreateThreadSendIndices(int threadId);
      int CreateThreadRecvIndices(int threadId);
      int UnpackSimpleRecvBuffers(int threadId);
      int PackSimpleSendBuffers(const state::base &inX, int threadId);
      void SetRemoteHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents)
      { 
        remoteHaloExtents = haloExtents;
        if(!localPartitionExtent.empty()){
          int numDim = haloExtents.size()/2;
          remoteHaloBufferExtents.resize(2*numDim);
          for(int iDim = 0;iDim < numDim;iDim++){
            pcpp::IndexIntervalType &leftHaloExtent(remoteHaloExtents[2*iDim]);
            if(leftHaloExtent.NNodes() > 0){ // Left remote halo for iDim direction
              int haloSize = leftHaloExtent[iDim].second - leftHaloExtent[iDim].first + 1;
              pcpp::IndexIntervalType remoteHaloBufferExtent(localPartitionExtent);
              remoteHaloBufferExtent[iDim].first = 0;
              remoteHaloBufferExtent[iDim].second = haloSize - 1;
              remoteHaloBufferExtents[iDim*2].Copy(remoteHaloBufferExtent);
            }
            pcpp::IndexIntervalType &rightHaloExtent(remoteHaloExtents[2*iDim+1]);
            if(rightHaloExtent.NNodes() > 0){ // Right remote halo for iDim direction
              int haloSize = rightHaloExtent[iDim].second - rightHaloExtent[iDim].first + 1;
              pcpp::IndexIntervalType remoteHaloBufferExtent(localPartitionExtent);
              remoteHaloBufferExtent[iDim].second = localBufferSizes[iDim] - 1;
              remoteHaloBufferExtent[iDim].first  = remoteHaloBufferExtent[iDim].second - haloSize + 1;
              remoteHaloBufferExtents[iDim*2+1].Copy(remoteHaloBufferExtent);
            }
          }
        }
        //         if(!localPartitionExtent.empty()){
        //         }
      };

      void SetThreadExtent(int myThreadId,pcpp::IndexIntervalType threadInterval)
      { 
        threadExtents[myThreadId] = threadInterval;
        if(!localPartitionExtent.empty()){
          int numDim = threadInterval.size();
          pcpp::IndexIntervalType &threadBufferExtent(threadBufferExtents[myThreadId]);
          threadBufferExtent = (threadInterval - globalPartitionExtent);
          threadBufferExtent += localPartitionExtent;
          for(int iDim = 0;iDim < numDim;iDim++){
            if(threadBufferExtent[iDim].first == localPartitionExtent[iDim].first){
              threadBufferExtent[iDim].first  = 0;
            }  
            if(threadBufferExtent[iDim].second == localPartitionExtent[iDim].second)
              threadBufferExtent[iDim].second = localBufferSizes[iDim] - 1;
          }
        }
      };
      void SetNumThreads(int numThreadsIn) 
      { 
        numThreads = numThreadsIn; 
        threadExtents.resize(numThreads,globalPartitionExtent);
        threadBufferExtents.resize(numThreads,localPartitionExtent);
        threadSendIndices.resize(numThreads);
        threadSendBufferIndices.resize(numThreads);
        threadRecvIndices.resize(numThreads);
        threadRecvBufferIndices.resize(numThreads);
        threadHaloBufferIndices.resize(numThreads);
      };
      int NumThreads() { return (numThreads); }; 
      // ----------------------------------------
  
      void SetLocalHaloExtents(const std::vector<pcpp::IndexIntervalType> haloExtents)
      { 
        localHaloExtents = haloExtents;
        if(!localPartitionExtent.empty()){
          // assert(!bufferSizes.empty());
          int numDim = haloExtents.size();
          localHaloBufferExtents.resize(numDim);
          for(int iDim = 0;iDim < numDim;iDim++){
            pcpp::IndexIntervalType &threadSendInterval(localHaloExtents[iDim]);
            // 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]);
            localHaloBufferExtents[iDim].Copy(subThreadSendExtent);
          }
        }
      };

      void CreateHaloBuffers();
      int CreateSendBuffers(const std::vector<pcpp::RemoteCollisionType> &inCollisions); 
      int CreateRecvBuffers(const std::vector<pcpp::RemoteCollisionType> &inCollisions); 

      int CreateSimpleSendBuffers();
      int CreateSimpleSendIndices();
      int CreateSimpleRecvBuffers();
      int CreateSimpleRecvIndices();

      int CreateMessageBuffers(int numComponents);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const double *sourceBuffer);
      int PackMessageBuffers(const int messageId,const int componentId,
                             const int numComponents,const double *sourceBuffer,
                             const int threadId);
      int SendMessage(const int messageId,const std::vector<int> &neighborRanks,
                      CommunicatorType &inComm);
      int ReceiveMessage(const int messageId,const std::vector<int> &neighborRanks,
                         CommunicatorType &inComm);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,double *targetBuffer);
      int UnPackMessageBuffers(const int messageId,const int componentId,
                               const int numComponents,double *targetBuffer,
                               const int threadId);
        
      void DestroySendBuffers();
      void DestroyHaloBuffers();    
      void DestroyRecvBuffers();
      int UnpackReceiveBuffers();
      int UnpackSimpleRecvBuffers();
      int PackSendBuffers(const state::base &inState);
      int PackSimpleSendBuffers(const state::base &inState);
      int PackSimpleSendBuffers(double inA,const state::base &inX,const state::base &inY);
      int Send(CommunicatorType &inComm);
      int SimpleSend(std::vector<int> &neighborRanks,CommunicatorType &inComm);
      int PostReceives(CommunicatorType &inComm);
      int PostSimpleReceives(std::vector<int> &sourceRanks,CommunicatorType &inComm);
      int CompleteReceives(CommunicatorType &inComm);
      int CompleteSimpleReceives(CommunicatorType &inComm);
      int Wait(CommunicatorType &inComm);
      void SetNeighbors(const std::vector<bool> &inNeighbors){
        haveNeighbor = inNeighbors;
      };
      inline std::vector<double *> &RecvBuffers()
      { return(recvBuffers); };  
      inline std::vector<std::vector<double *> > &HaloBuffers()
      { return(haloBuffers); };
      inline std::vector<double *> &SendBuffers()
      { return(sendBuffers); };
      inline std::vector<std::pair<int,std::vector<size_t> > > &HaloRecvIndices()
      { return(haloRecvIndices); };
      inline std::vector<pcpp::RemoteCollisionType> &SendCollisions()
      { return(sendCollisions); };
      inline std::vector<pcpp::RemoteCollisionType> &RecvCollisions()
      { return(recvCollisions); };
      inline double *HaloBuffer(int fieldIndex,int directionIndex)
      { return(haloBuffers[fieldIndex][directionIndex]); };
      inline std::vector<pcpp::IndexIntervalType> &ThreadExtents() { return(threadExtents); };
      inline std::vector<int> &StateFieldIndices() { return(stateFieldIndices); };
      inline pcpp::IndexIntervalType &GlobalExtent() { return(globalGridExtent); };
      inline pcpp::IndexIntervalType &PartitionExtent() { return (globalPartitionExtent); };
      inline std::vector<pcpp::IndexIntervalType> &RemoteHaloExtents() { return(remoteHaloExtents); };
      inline std::vector<pcpp::IndexIntervalType> &LocalHaloExtents() { return(localHaloExtents); };
      const std::vector<pcpp::IndexIntervalType> &RemoteHaloBufferExtents() const { return(remoteHaloBufferExtents); };
      const std::vector<pcpp::IndexIntervalType> &LocalHaloBufferExtents() const { return(localHaloBufferExtents); };

      std::vector<messagebuffers> &CommunicationBuffers() { return(communicationBuffers); };
      const std::vector<std::vector<size_t> > &SendIndices() const {return(sendIndices);};
      const std::vector<std::vector<size_t> > &RecvIndices() const {return(recvIndices);};
      inline std::vector<int> &PeriodicDirs() { return(periodicDirs); };

      void ReportSimpleBufferContents(std::ostream &outStream);
      ~halo(){
        DestroyAll();
      }
      void DestroyAll()
      {
        DestroyHaloBuffers();
        DestroySendBuffers();
        DestroyRecvBuffers();
      };
    protected:
      
      bool ownSendBuffers;
      bool ownRecvBuffers;
      bool ownHaloBuffers;
      bool haveRecvData;
      bool haveSendData;
      bool doFill;
      int numThreads;
      size_t numPointsBuffer;
      //    int numStateComponents;

      std::vector<int> periodicDirs;

      // These are the extents of the partition and
      // halos, resp.
      pcpp::IndexIntervalType globalGridExtent;
      pcpp::IndexIntervalType globalPartitionExtent;
      pcpp::IndexIntervalType localPartitionExtent;
      std::vector<size_t> localBufferSizes;

      std::vector<bool> haveNeighbor;
      std::vector<pcpp::IndexIntervalType> remoteHaloExtents;
      std::vector<pcpp::IndexIntervalType> localHaloExtents;
      std::vector<pcpp::IndexIntervalType> localHaloBufferExtents;
      std::vector<pcpp::IndexIntervalType> remoteHaloBufferExtents;
      std::vector<pcpp::IndexIntervalType> sendExtents;

      std::vector<pcpp::IndexIntervalType> threadExtents;
      std::vector<pcpp::IndexIntervalType> threadBufferExtents;

      // For each send required, these indicate the index ranges
      // and a flat array of indices into the state data
      std::vector<pcpp::RemoteCollisionType> sendCollisions;
      std::vector<std::vector<size_t> > sendIndices;
      std::vector<std::vector<std::vector<size_t> > > threadSendIndices;
      std::vector<std::vector<std::vector<size_t> > > threadSendBufferIndices;
      
      // These indicate the receive buffer intervals
      std::vector<pcpp::RemoteCollisionType>            recvCollisions;
      std::vector<std::vector<size_t> >                 recvIndices;
      std::vector<std::pair<int,std::vector<size_t> > > haloRecvIndices;
      // receive intervals for threads
      std::vector<std::vector<std::vector<size_t> > >   threadRecvIndices;
      std::vector<std::vector<std::vector<size_t> > >   threadRecvBufferIndices;
      std::vector<std::vector<std::vector<size_t> > >   threadHaloBufferIndices;


      
      
      // Communication buffers are densely stuffed and
      // packed/unpacked with special functions
      std::vector<messagebuffers> communicationBuffers;
      std::vector<double *> sendBuffers;
      std::vector<double *> recvBuffers;
      
      std::vector< std::vector<double *> >  haloBuffers;
      
    };
    
    template<typename StateType>
    void ConvertLegacyState(int numDim,size_t numPoints,StateType &inState)
    {

      // convert density
      int rhoIndex = inState.GetFieldIndex("rho");
      if(rhoIndex >= 0){ // then rho already exists, just switch its buffer
        double *legacyDataBuffer = inState.template GetFieldBuffer<double>("cv01");
        inState.SetFieldBuffer("rho",legacyDataBuffer);
      } else { // otherwise rename cv01 to rho
        inState.RenameField("cv01","rho");
      }
      
      // convert energy
      std::ostringstream Ostr;
      Ostr << "cv0" << numDim+2;
      const std::string legacyRhoE(Ostr.str());
      int rhoEIndex = inState.GetFieldIndex("rhoE");
      if(rhoEIndex >= 0) { // rhoE already exists, buffer switch
        double *legacyDataBuffer = inState.template GetFieldBuffer<double>(legacyRhoE);
        inState.SetFieldBuffer("rhoE",legacyDataBuffer);
      } else {
        inState.RenameField(legacyRhoE,"rhoE");
      }

      // convert momentum
      int rhoVIndex = inState.GetFieldIndex("rhoV");
      if(rhoVIndex < 0) {
        inState.AddField("rhoV",'n',numDim,8,"momentum");
        inState.Create();
      }
      double *targetData = inState.template GetFieldBuffer<double>("rhoV");
      for(int iDim = 0;iDim < numDim;iDim++){
        std::ostringstream Ostr2;
        Ostr2 << "cv0" << iDim+2;
        const std::string fieldName(Ostr2.str());
        double *srcData = inState.template GetFieldBuffer<double>(Ostr2.str());
        std::memcpy(&targetData[iDim*numPoints],srcData,sizeof(double)*numPoints);
      }

      // convert auxVars
      const std::string pc2Name("scalarVars");
      int scalarIndex = inState.GetFieldIndex(pc2Name);
      
      int numComponents = 0;
      const pcpp::field::metadataset &legacyMeta(inState.Meta());
      pcpp::field::metadataset::const_iterator metaIt = legacyMeta.begin();
      while(metaIt != legacyMeta.end()){
        const pcpp::field::metadata &fieldMeta(*metaIt++);
        const std::string &fieldName(fieldMeta.name);
        std::string::size_type x = fieldName.find("aux");
        if(x == 0){ // this is an aux var
          numComponents++;
        }
      }
      if(numComponents > 0){

        if(scalarIndex < 0){
          inState.AddField(pc2Name,'n',numComponents,8,"massfraction");
          inState.Create();
        }
      
        targetData = inState.template GetFieldBuffer<double>(pc2Name);
        for(int iComponent = 0;iComponent < numComponents;iComponent++){
          std::ostringstream nameStr;
          nameStr << "aux" << (iComponent < 9 ? "0" : "")
                  << iComponent+1;
          const std::string legacyName(nameStr.str());
          double *srcData = inState.template GetFieldBuffer<double>(legacyName);
          size_t targetOffset = iComponent*numPoints;
          std::memcpy(&targetData[targetOffset],srcData,sizeof(double)*numPoints);
        }
      }
    };

  }; 
};
#endif